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SUMMARY 


The operation of the MAST code, which computes transient solutions to the 
multiphase flow equations applicable to all-speed flows, is described. Two-phase flows are 
formulated based on the Eulerian-Lagrangian scheme in which the continuous phase is 
described by the Navier-Stokes equation (or Reynolds equations for turbulent flows). 
Dispersed phase is formulated by a Lagrangian tracking scheme. The numerical solution 
algorithms utilized for fluid flows is a newly developed Pressure- Implicit algorithm based 
on the operator-splitting technique in generalized non-orthogonal coordinates. This operator 
split allows separate operations on each of the variable fields to handle pressure-velocity 
coupling. The pressure correction equation thus obtained has the hyperbolic nature and is 
effective for Mach numbers ranging form the incompressible limit to supersonic flow 
regimes. The present code adopts a non-staggered grid arrangement, thus the velocity 
components and other dependent variables are collocated at the same grid. A sequence of 
benchmark-quality problems, including incompressible, subsonic, transonic, supersonic, 
gas-droplet two-phase flows, as well as spray-combustion problems, were performed to 
demonstrate the robustness and accuracy of the present code. 


CHAPTER 1 


INTRODUCTION 

The partial differential equations governing fluid flows in a typical liquid-fueled 
rocket engine are nonlinear and are too complex to obtain the solution in analytical form in 
general, therefore computational methods for simulating flow fields must be employed. 
Since problems of practical interest include a very wide range of flow situations, it would 
be desirable to develop a single numerical algorithm capable of reliably handling many flow 
conditions in an efficient manner. 

The numerical modeling of the physical mechanisms of turbulence effects, 
multiphase flows, and combustion are very complex individually and the formulation of an 
algorithm combining all of the above mechanisms would be a formidable task. However, 
an approach can be adopted such that one aspect of the numerical simulation can be 
considered individually and collectively common to the above process; this aspect being the 
transient all-speed Navier-Stokes flow solver. As the first logical step in attaining a fully 
comprehensive chemical reacting flow simulator, we develop a versatile and generalized 
flow solver which could efficiently handle multidimensional, incompressible and 
compressible, laminar and turbulent, reacting and non-reacting, steady and unsteady, 
single- and multi-phase flows. This report is intended as reference guide for users of the 
MAST (Multiphase All-Speed Transient) family codes. 

The primary features of the MAST code can be summarized as follows: 

1. A time-accurate formulation is employed, thus a steady state solution is achieved 
through a temporal marching approach. 

2. A pressure-based method using pressure and velocity as primary variables is 
formulated on a non-staggered, general, non-orthogonal coordinate system. 

3. The pressure-velocity coupling is handled through a modified, non-iterative 
operator- splitting algorithm valid for all speed flow situations. 
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4. The set of discretized equations is solved by using a conjugate gradient (CGS) 
method. 

5. The dispersed phase is handled through a Lagrangian tracking technique 
embedded in the pressure-velocity predictor-corrector algorithm for efficient two- 
way coupling. 

In the following, fundamental elements for establishing the MAST code will be 
described. These are followed by a series of bench-mark fluid dynamics calculations. A 
user's guide is also provided in Appendixes. 
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CHAPTER 2 


GOVERNING EQUATIONS 

The unsteady, Newtonian, viscous compressible flow is governed by the following 
Navier-Stokes equations: 
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In the above p is the density, Uj the velocity components, p the pressure. Si are body 
forces and are the components of deviatory stress tensor. In equation (3), h is the 
stagnation enthalpy of the fluid, and is related to the specific energy e (which is a function 
of temperature T) by 

h = e(T) + — +— uni; 
p 2 

The viscous constitutive relation 
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and the density distributions are provided by the equation of state. 
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P 
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To treat arbitrary boundaries, these governing equations are transformed to a 
general form based on a nonorthogonal coordinate system (£,T|,Q. There are three ways to 
represent velocity field on a nonorthogonal grid system; i.e. (1) Cartesian Velocity 
Components [1,2], (2) Contravariant velocity component [3] and (3) resolute or velocity 
projections along the coordinate directions [4,5]. The present approach uses the Cartesian 
components of the velocity vector stored at grid points. This has the advantage of excluding 
the so called extra "curvature" terms in the momentum equations. Using the standard 
transformations formula [6], the continuity equation becomes 

f4[| (pu)+ i (pv)+ 4H =0 (6) 

and the remaining governing equations for the dependent variables b in the general 
coordinate system can be expressed in a compact form as follows: 

^(p4>)+y ^(P u <f)+^(P v '»)+4(p w 4.) 

^ / 2 2 2 \ ^ 3 4/ 1 2 2 \ 

~ 9^[t( 8ii+ gl2+ g13 ) a^J + aii |r( g21+ g22+ g23 )anJ + 

y-(g31+g32+g33)|^ +S<1> ( 7 ) 



and and S'** are the associated diffusivity and source terms for the variable <)> ( = u, v, w, 
T, k, e, etc. ). The detailed expressions of the source terms and the cross derivative terms 
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due to grid non-orthogonality are listed in Table 1. The geometric coefficient gij and the 
Jacobian J are defined as 


g ij “ 
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For turbulent flow computations, the MAST code has employed the two-equation 
K-e model with wall-functions. In this case, the molecular viscosity and thermal 
diffusivities are augmented by the effective viscosity and effective diffusivities as: 


M-eff = + 
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The transport equations for the k and e can also be cast in the standard form with 
source terms listed in Table 2. 

COMPUTATIONAL METHODS 

A. Discretization 


Discretization of the governing, differential equations uses the finite volume 
approach [7]. Differencing in the temporal domain employs the Implicit Euler scheme. All 
the dependent and independent variables are stored at the same grid location and variables at 
the finite control volume boundaries are interpolated between adjacent grid points. In the 
spatial domain, the diffusion terms are approximated by the central differencing scheme. To 
handle sharp shock-capturing without oscillatory fields near the shock discontinuity, a high 
accuracy scheme for handling the convection terms is implemented. This scheme is briefly 
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described for the convective flux term in the x direction. Let the discretized formulation of 
the invicid component for grid point be 


afruo), af, f 4~ f 4 

3x 1 3x 1 Ax 

then the control- volume interface flux is computed according to: 
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In the above, h ( represents a first - order numerical flux, and can be computed from: 
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The flux-limited values of df are computed as follows: 


df + i4 = p. l u + i((t>i +I -<t>i) 
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In contrast to the minimod operator used in [8] , we use the following "compression" 
parameters determined by a local normalized pressure gradient monitor as follows: 
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where p i = min(gi , g i+1 ) 
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where 


0 


Si = 


(A + Boij) 

Acq +B 
A + B 

A + Baj 


; a, < 0 

; 0 < (Xi < 1 


; a* >1 


and aj = ^±1 — — 
Pi - Pi-i 


In this study, A and B are set to be 1 and O is chosen to have high resolution for specific 
flow regimes. Suggested values are given in the results sections. 


B. Velocity-Pressure Coupling 


A typical computational procedure for compressible flows employs density as a 
primary variable and extracts pressure from an equation of state. This density-based 
method has difficulties in handling incompressible or low-Mach number flows, due to 
weak pressure-density coupling. 

Computation for incompressible flow is generally performed using a pressure based 
method which is characterized by the use of pressure as one of the primary dependent 
variables. However, the absence of an explicit equation for pressure imposes that the 
pressure can only influence the velocity field through the momentum conservation with the 
continuity as a compatibility condition. Hence provision of algorithm for the velocity- 
pressure coupling is essential. There are several algorithms employing a semi-implicit 
iterative coupling procedure; e.g. SIMPLE [7] and its variants [9]. In the MAST code, an 
efficient non-iterative method based on the operator- splitting technique is extended to non- 
orthogonal grid coordinate. The discretized governing equations are solved in a time- 
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marching fashion. In each time step, a fixed number of correction stages follow an initial 
predictor stage, and this predictor-multiple corrector procedure yields a time-accurate result 
As is the case with the basic operator-splitting method depicted in [10], we seek to 
discretize the governing equations as the following by using the simple Euler implicit 


scheme in temporal domain: 

i{(pui) n+1 -(pui)"} = H‘(uj n+1 )— Aj p n+1 + Si (19) 

^{(ph) n+1 -(ph) n ] = J'(h n+1 ) + ~ ^g~ ~ + Qi (20) 

In Eq. (19) and (20), finite difference operator H'( ) and J’ ( ) contains the convection and 
diffusion operations. They can be further split into: 

H’(ui n+1 ) = A 0 ui( )+H(u/ )) (21.a) 

J'(h n+1 ) = B 0 h ( K j(h ( ^ (21.b) 


in which A 0 and Bo denote the finite difference coefficient matrix for diagonal elements. 
The splitting of operator in (21. a) and (21. b) allows different intermediate time steps be 
used. This technique was introduced to derive the velocity-pressure coupling sequence in 
several intermediate steps. Let superscripts *, ** and *** represent predictor, first and 
second corrector field values obtained during the splitting procedure. For simplicity, only 
single-phase, laminar pressure-velocity coupling equations are described here first 
(a) Predictor step. The current velocity, density and pressure field are used to calculate 
finite-volume discretization coefficient to solve the implicit momentum predictor step: 


f _n _n-l n 

— -A 0 n* =H(ui)-Aip n +Si+ ■ — 


St 


In this step, H'(u *j) = A 0 u*j +H(u* i ) 


( 22 ) 


(b) First momenmm corrector step. To derive a pressure-velocity coupling formulation 
satisfying conservation of mass is the goal of this step. The form of pressure equation 
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depends greatly on the discretization of the continuity equation. In this study, we propose 
the following form of the continuity equation: 


p*— p n 
8t 


Ai(p n Uj* *) + A i (p'u i *) = 0 


(23) 


where p'=p*-p n 


This equation differs from the original PISO formulation and is consistent with the 
linearization procedure used in momentum equations as well as other transport equations 
To utilize this equation, the momentum equation (2) can be discretized in an explicit form 
for the second corrector form as: 


f 





Uj* * 


H’(uj*)- Aj p* + Si+- 


>" +1 Uj" 
6t 


(24) 


Note that in this step, H'(uj * *) = A 0 u i ** + H(u i *) 

We now subtract Eq.(22) from Eq. (24) and then take the divergence of the resulting 
equation and substitute into Eq. (23), which gives the implicit pressure correction equation 
for the first corrector step: 


,8tRT 


+ A i - Ai (P" D i 4 i )}(P* -P" ) = - A i ( pD u i*) 


(25) 


where 


Dj = 


V 

J7 -a ° 


Y-l 

) 


Note that this equation involves a convection term which enhances the convergence and 
properly takes into account the hyperbolic nature of supersonic flows, and enables 
capturing shock waves. In Eq. (22), p* has been substituted by the equation of state 
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The energy field for the first correction step is then solved as: 


8t 


-B 0 


o" 1 n * ^ 

h* = j(h*)+£^-h"+£-^-+Q i * 


( 26 ) 


(c) Second corrector step. The first corrector step satisfies the continuity equation (23). To 
satisfy momentum and other equations further, a second pressure correction equation was 
derived using an explicit second momentum corrector and a discretized continuity equation: 
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where p”=p**-p* 

The resulting form of the second pressure correction equation becomes: 





+ Aj 



RT*y 


u, 


** 


(29) 


** 

IT ** p 
Here, p = 

RT 

was used for deriving Eq. (29). This second pressure-correction equation also differs from 
the PISO algorithm by a convection term on the left hand side and an extra term on the right 
hand side. The energy equations are solved again to obtain h** and eventually T** . It is 
shown in [1 1] that the errors introduced by the splitting procedure in the temporal domain 
vanish with 8t and do so at least as fast as the temporal discretization enror introduced by 
the Euler discretization scheme. The variables solved after the second corrector steps are 
then taken to be the field values at the next time level, n+1. It should be noted that, for 
incompressible flows for which p is a constant, the present method reduces to the original 
PISO algorithm [10]. 
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C. The Collocated Grid Arrangement 


It is known that the non-staggered grid system can produce undesirable checker- 
board (oscillatory) phenomena in the pressure field, due to the decoupling effect between 
velocities and pressures. As suggested by Rhie [1], an explicit fourth order pressure 
damping term may be added to the right hand side of pressure increment equations to damp 
out pressure oscillations. This approach can be further modified to include the transient 
term for a time-marching formulation. However, the steady state solution obtained in this 
way is found to be strongly influenced by the size of time step chosen, which is 
undesirable [26]. 


The pressure-velocity coupling procedure described above is general and is not related to 
any particular grid arrangements and spatial discretization schemes. In this study, this 
coupling algorithm has been formulated in a general curvilinear non-orthogonal body-fitted 
coordinate system with a non-staggered grid arrangement, i.e. ,all variables are located at 
the same grid point. A new weighted intetpolation method was developed to avoid checker- 
board pressure field frequently encountered in pressure-based method utilizing non- 
staggered grid. This treatment is briefly described here. More detailed justifications can be 
found in [11]. For simplicity, let the finite difference operator acting on the u velocity 
component (in the x direction) at grid points i and i+1 be written as: 


A i u i =H i -|£| i +S 


; , (P U )i n 

' St 


(30) 


A i-H u i+1 ~ H i+1~ |^|i+l + S i+1 + ~' g t ~~ @1) 

In the control-volume based finite difference method used in this study, the main concern is 
the estimation of control-volume face value u. In the adaptive interpolation method 
developed in [1 1], the face value is evaluated according to: 


u =p t u 
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where ft x =l-max(g i+1 , &) 
1 +- 
2 



in which g* = 


Pi+i+Pi-i-2pi 

Pi+i+Pi-i+^Pi 


In Eq. (32), u i+i/2 is simply an arithmetic averaging , =(ui +ui+i )/2, and u i+i /2 is 
estimated by the following approximate operator: 


(Ai+A i+1 )u j =(A i+1 -A i )(u i+1 -u i ) + H i +H i+1 

in — 



l + Sj+ Sj + i+ 

1 +- 


2(P u ) n i 

i+— 

2 

8t 


(33) 


Detailed testings of this method and the comparisons with other non-staggered treatments 
[12,13] can be found in [11]. 


D. Equation Solver 

The system of finite difference approximation equations formed by equations (22), 
(25), (26) and (29) produces a nine-banded matrix. To enhance the matrix solver 
efficiency, the conjugate-gradient-squaied (CGS) method as described in [14] is employed. 

DISCRETE PARTICLE TRACKING METHOD 


In the MAST code, the dispersed phase, i.e. the spray, is considered to be 
composed of discrete computational parcels, each of which represents a group of droplets 
of similar size, velocity, temperature etc. The distribution function in droplet size, velocity, 
temperature, spray pattern produced by the fuel injector (or atomizer) is statistically 
sampled and the resulting computational parcels are followed as they interact and exchange 
moment, energy and mass with the surrounding gas. The present numerical model employs 
this Monte Carlo method and the basic ideas of implementing this technique in the present 
CFD method are presented in Ref [15], which is reproduced in Appendix C. Detailed 
descriptions of the sampling technique of specified distribution functions for characterizing 
the initial conditions for spray and the subsequent droplet-fluid interactions can also be 
found in Ref [15] and [16]. 
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One of the important aspects in spray combustion modeling is the dense spray 
effects which include atomization process, drop breakup, droplet collision and coalescence. 
Atomization process occurs on time and length scale too short to be resolves with practical 
computation grid sizes and time steps. Thus, atomization should be modeled as a sub-grid- 
scale process. To account for the dense spray effects, the present study employs the drop 
collision and coalescence model [33] and the Taylor analogy breakup (TAB) model [34]. In 
the drop collision model, the probability distributions governing the number and outcomes 
of the collisions between two drops are sampled randomly in consistency with the 
stochastic particle tracking method. The TAB model utilizes an analogy between an 
oscillating and distorting droplet and a spring-mass system. The present breakup model is 
based on the reasonable assumption that atomization and drop breakup are indistinguishable 
processes within a dense spray near the nozzle exit. Accordingly, atomization is prescribed 
by injecting drops which have a characteristic size equal to the nozzle exit diameter. 


CHEMICAL REACTIONS 

In the current version of the MAST code, it is assumed that the spray is dilute and 
the liquid fuel droplets act as distributed sources of fuel which evaporate to form a cloud of 
vapor. The combustion process thus is treated as gaseous diffusion flames. Here, an 
idealized approach for physically-controlled diffusion flames is to invoke a fast-chemistry 
assumption which the chemistry is sufficiently fast and intermediate species do not play a 
significant role. 

In cases where chemical reactions are assumed to be in equilibrium, the equilibrium 
energy minimization procedures, the reader can refer to the well-known NASA Lewis 
CEC76 program [17]. A simplified version of equilibrium chemistry procedure based on 
[17] is incorporated into the MAST code for calculating mixture compositions, 
thermodynamics and transport properties (also refer to JANAF table [18]). 

In the turbulent diffusion flame model, the influence of turbulence on combustion is 
taken into account by relating the fluctuations of mass fractions. This implies that fuel and 
oxidizer can coexist in the same place but at a different time. The effect of turbulent eddies 
on thermodynamical properties is modeled via the introduction of the probability function 
P(£»Xi) [20]. This function contains information of both mean (f) and variance of the 
mixture fraction [(f-f ) 2 ]=g. There variable f and g can be obtained by solving the transport 
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equations. The density- weighted mean value (^) of any property are evaluated by 
convolving the property functions with the pdf, P(£,xi): 

<i> = /<l>(S)P(txi)dS 

Numerous pdfs are available in the literature. The MAST code adopts the (3-pdf which is 
known as the widely applicable one [21], 
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Chapter 3 


THE COMPUTER PROGRAM 

A. General Structure 

The MAST family computer program consists of a set of subroutines controlled by 
a short main program. The fundamental structure is illustrated in Figure 3.1, showing a 
top to bottom block diagram encompassing the entire calculation procedure. Table 3-1 gives 
a brief description of the subroutines found in the MAST code. List of important 
FORTRAN symbols are given in Appendix A. In the following , the input file and the 
handling of boundary conditions are described. 

The Input structure used in the MAST code utilizes blocks of optional namelists 
under several keywords. Each of the keywords and associated namelists are described 
here. These "keywords" are optional, i.e., skip if not needed. 

There are ten keywords built into the current MAST code. These are : GRID, 
BOUND, SOLV, PROPERTY, TURBULEN, SPRAY, REACTION, RUN and 
ENDJOB. In addition, there is one extra block called CONTROL Block used for 
identifying numerics options in the code. The CONTROL block is usually created first in 
the input file to be read into the program through Login Unit 1. All entries are optional. If 
certain entries require numerical values or logical values, they are entered after the entry 
names separated by at least one blank space. Block name, entry names, the possible range 
of values and a short description of the variable are described as follow. 

1 . CONTROL Block : This block is always put at the top of the input file. 

RESTART: It activates the usage of restart file as initial 

conditions by reading through LU=4. 

SWIRL: It activates the swirl velocity calculations. 

IMON, 

JMON: It specifies the monitor point at (IMONJMON), 

The default value is (2,2). 

Only one variable can be specified as 
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MONU,MONV, 


MONTEMP, 

MONTK, 

MONTE: 

ERRCG: 

NCGM: 

ERRM: 

INCOMP: 

COMPRES: 

OMGD: 

NCRT: 

OMGF: 


OMGT: 

PHI: 

OMGPHI: 

END: 

2. GRTD Block : 

NX: 

NY: 

XLEN: 

YLEN: 

READXY: 


a monitor variable tracking on screen. 

Such as velocity at (x,y) location, pressure, 
temperature, turbulent kinetic energy, 
dissipation rate. The default is MONU. 

Convergency criterion for conjugate gradient 
Matrix solver. The default value is 0.01. 

Maximum CG solver iteration number. 

The default number is 100. 

Termination criterion for steady state solutions 

using the time marching scheme. The default value is 0.0001 . 

Incompressible flow calculations (=1). 

(INCOMP=0) Compressible flow options. 

1 for supersonic flows, 0 (default) for low speed 
flows. 

Number of corrector steps (default=2). 

Interpolation parameter for face velocities 
1 for interpolation and 0 for averaging. 

The default value is 1. 

For supersonic temperature field relaxation. 

Parameter in the finite difference limiter. 

Parameter for Upwinding scheme. 

End of Block input. 


Number of grid points in the x direction (>3). 

Number of grid points in the y-direction (>3). 

It identifies the length of the calculation domain in the x direction. 
The length of the calculation domain in y direction. 

It activated reading of grid point X(I,J), 
and Y(I,J) from input file "sgrid.d". The input 
format is given as: 

READ(4,9910) NX, NY 
READ(4,9920) ((X(I,J),I=1,NX),J=1,NY) 
READ(4,9920) (( Y (I, J),I= 1 ,NX), J= 1 ,NY) 

9910 FORMAT(2I5) 
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UNIFORM: 

9920 FORMAT(10E10.4) 

It specified a uniform grid system. 

XDIR,YDIR: 

It specifies an X-direction grid (only 

1ST: 

dependent on I) or a Y-direction grid (only dependent 
on J).Either one must be put at first of a line. 

The grid cell starts from 1ST index 

IEND: 

The grid cell ends with IEND index 

DST: 

The grid cell starts from DST 

DEND: 

The grid cell ends with DEND 

EXP: 

stretching factor 

DELT: 

first grid cell size. EXP and DELT can 


be specified only once. 


"power law" * Such a grid generation can be described as follows: 

(a) Explicit "power law": 

For example, XDIR 1ST 10 IEND 25 DST 3.2 DEND 5.6 EXP 1.5 The power-law formulation for 
such a grid distribution is, 

( t tct \ EXP 

X(I, J) = DST+ ■ — (DEND- DST) 

VlEND- ISTj 


I=IST, IEND for EXP>0 

X(I,J) = DEND- 


IEND-I V™ 


IEND- 1ST ) 


(DEND- DST) 


I=IST IEND for EXP<0 

. If EXP = ±1, grid is uniform 

. If EXP > 1, grid is compressed close to DST and expanded close to DEND. 

. If EXP < -1, grid is expanded close to DST and compressed close to DEND. 
. EXP >1 and exp < - 1 can be used to give a symmetric grid distribution. 


(b) Implicit "power-law" : 

If the user specified first grid cell size DELT instead of EXP, the grid generator 
automatically computes the stretching factor, EXP, and distributes the grid maintaining the 
first grid size DELT = X(IST +1, J) -X(IST, J) for DELT > 0 or the last grid size IDELTI 
=X(IEND, J)-X(IEND-1, J) for DELT < 0 


3. Property block (PROPi 
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VISCOS: 

It specifies the fluid viscosity. The default 
value is 1. 

DENGAS: 

It specifies the fluid density. The default 
value is 1. 

CPGAS: 

It specifies the fluid specific heat The 
default value is 1. 

KGAS: 

It specifies the fluid thermal conductivity. 
The default value is 1. 

PRG: 

It specifies the fluid Prandtl number. The 
default value is 0.74. 

TEMP: 

It specifies the flow field initial temperature. 
The default value is 0. 

OMEGA: 

SMW 

Angular momentum 

GAMMA 

= Cp 

^ C 
'-v 

K: 

conductivity of the Block region 

DENSITY: 

density 

GEN: 

Heat generation 

UIN: 

Initial field of U 

4.Runtime block (RUN'! 

DT: 

It specifies time step. 

TIME: 

Maximum time to be computed. DT and TIME can 
be specified only once. 

TIME =DT *NTIME 

NTIME: 

Maximum time step number to be computed. 

NPR1: 

It specifies the screen monitor output frequency. 

NEX: 

It specifies the example number. The 
user can code in SUBROUTINE EXAMPL. 

VIN: 

Initial field of V 

STOP: 

If this keyword is specified, MAST code 

only checks the input data and no progress is going 
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U,V: 

Velocity field will be solved. 

P 

Pressure will be solved. 

TEMP: 

Temperature will be solved. 

TK, TE: 

Turbulence kinetic energy and dissipation rate 
will be solved. 

SW: 

Swirl velocity will be solved. (SWIRL must be 
activated in Selection block). 

PATC: 

Particle tracking is active. 
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1ST, END, If 1ST = END or JST = END, it 

JST, END : becomes a line boundary condition. 

Boundary conditions specified always take such a form : 

1ST, END, JST, END, B.C. type, variable. 

The following keywords used as B.C. type : 

INLET: Inlet boundary condition, followed by variables 

OUTLET: Outlet boundary condition 

WALL: No split wall boundary condition, followed by 

variables. Wall function of turbulence model activated. 
SYMMETRY : It specifies a symmetric b.c. 

CYCLE : It specifies a cyclic b.c. 

BLOCK : It specifies a blockage. Wall function 

of turbulence model activated on block face. 


Variable keywords (specify b.c.O) 


U : x; -direction velocity 

V : y-direction velocity 

TK : Turbulence kinetic energy in INLET 

TE : Turbulence dissipation rate in INLET 

TEMP :: Temperature in INLET WALL or BLOCK 


9. Reaction block (REACTIONS 


NSPE 

NELE 

NFROS 

FROSEN 

HP 

TP 

ELEMCONS 


CHEMFORM 


: number of species to be solved. 

: number of elements. 

: number of species to be frozen if frozen is activated. 
: activate ffosen. 

: assigned enthalpy and pressure. 

: assigned temperature and pressure. 

: element conservation constants, repeat nele times. 

for example: air 

1 O 2 + 3.76 N 2 

number of element 0= 1 *2=2 

number of element H=3.76*2=7.52 

the conservation constants are 2,7.52, or 1,3.76 

: chemical formulation of species. 
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input format (4X, A8, 20F2.0), repeat NSPE times. 

A8 chemical formulation 

F2. 0 number of individual element presents in chemical formulation. 
10. ENDJOB : ENDJOB identifies the end of the input file. 
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CHAPTER 4 


TEST EXAMPLES 

In this chapter, various sample calculations are presented. These calculations 
include incompressible, compressible, invicid, viscous, laminar, turbulent flows. Two- 
phase flows and chemical reacting flows calculations are also included. 

1. Driven Cavity Flow 

A square 2-D cavity with the top wall moving at a constant speed is calculated. A 
non-uniform grid with finer grids near the wall and a Reynolds number of 1000 are used. 
The input decks written on the Logic Unit (LU) 1 are listed in Figure 4.1.1. This input file 
consists of the Selection Block and other six blocks. Any comments can be done following 
a semicolonO) at the beginning of a line or after name lists. For incompressible flow 
calculations, typical selections of parameters for interface interpolation and finite different 
schemes are given here. The U momentum was monitored at the grid point (4,3). The 
convergence criteria in the CG solver is 1. E-l and the calculations are terminated at 10-4. 

In the GRID Block, 51 x 51 grid was used. The grid was constructed such that the 
first I grid point coincides with x = 0, the 26th I grid point is at x = 0.5, and the 51st I grid 
point is at x = 1.0. The exponential expanding scheme was used for cluster finer grid near 
walls with rate 1.5 as illustrated in Figure 4.1.2. The same expansion ratio and grid 
arrangement are the same in the y-direction. In the BOUND block, boundary conditions 
were identified. The top wall was represented by I grid from 1 to 51 and J = 51, U is set to 
be 1. East and West side boundaries were set to no slip conditions (WALL). In the 
PROPERTY Block, VISCOS was set to 1.0E-3 to represent Re = 1000. As can be 
realized, the MAST code formulates the equation from a normalized set of variables. As 
long as characteristic variables used to normalize the governing equations are consistent, 
different combinations are possible. In the SOLV block, U, V, P were identified. In the 
RUN block, time step (DT) was set to 1.0 and 150 time steps NTIME were marched. 
NPR1 was set to 10 to activate printing monitoring results, ENDJOB gave the end of input 
file. 
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Part of the program output which is written to a monitor screen for this example is 
included here (see Figure 4.1.3). The convergence or evolution history is printed for each 
time step. The entry in the first column STEP is the current time step. The second entry, 
TIME is the normalized time. The entry FMON display the variable being monitored at a 
specified point. The fourth entry indicates the variable which has the greatest change from 
on time step to the next. The entry FERRM is the maximum rms change of any of the 
variables occurring in that time step, a measure of how well the steady state being 
approached. If this number reduced to 10~ 3 , the solutions essentially achieve steady state. 

Figure 4.1.4 shows the computed streamline contours. This case converges in one 
run in 83 CRAY/XMP seconds with no restarting necessary. The core memory used is 
0.39 M(64-bit) words. 

2. Two-D Circular Cylinder 

The flow over a two-dimensional circular cylinder is an example of the external 
unsteady flow when Reynolds number is not too small. A 41x41 O-grid which wraps 
around the cylinder is used for the grid. A simple algebraic grid generation was coded in 
Subroutine EXAMPLE, in this case NEX 2 was used to activate the example. The grid 
system is shown in Figure 4.2.1. A fixed velocity was set at the outer boundary and a 
cyclic boundary condition was set at the front center line. The input decks are listed in 
Figure 4.2.2. Stream function plots at several time steps are plotted in Figure 4.2.3. More 
detailed calculations of this case can be refereed to [22]. 500 steps of this calculation take 
175 CPU seconds and 0.29 M words core memory. 

3. Laminar Backward- Facing Step Flow 

The flow over a backward-facing step with a two-to-one expansion ratio was 
computed for a Reynolds number 800. A uniform 81x51 grid with blockage I from 1 to 
1 1 and J from 1 to 26 was used. The relevant parameters used in this calculation are given 
in Figure 4.3.1. An experimental and theoretical study object this problem has been 
provided by [15]. Stream function plots of the calculation is presented in Figure 4.3.2. The 
calculated reattachment length at the lower wall was about 11.7 step height. This result 
compared favorably with other second order accurate finite difference results. The 
experimental data of [15] is about 14 H for this case. This case took 300 CPU seconds on a 
CRAY/XMP for 350 time steps and requires 0.54 MW core memory . 
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4. Turbulent Backward Facing Step Flow 


A turbulent flow with Reynolds number 10,000 and three to two expansion ratio 
was calculated. A 59 x 35 non-uniform grid was used with an inlet mean velocity placed 
upstream of the expansion. The grid was clustered near all the solid walls using an 
exponentially stretching function. The standard high Reynolds number k-e model with wall 

functions was used in this calculation. The relevant parameters are given in Figure 4.4.1. 
The grids are shown in Figure 4.4.2. The calculated stream lines are shown in Figure 4.4.3 
and the calculated kinetic energy levels are shown in Figure 4.4.4. The predicted 
attachment length is 6.26 H. These results compared well with other results from well 
established codes. It took 60 CPU Seconds to run this case for 90 time steps and required 

0. 35 MW core memory.. 

5. Subsonic Channel Flow with a Circular Arc Bump. 

A series of two-dimensional invicid flow calculations following the same 
geometries used in the literature [23,24,25] have been chosen to demonstrate the invicid 
flow calculations. A circular arc bump thickness-to-chord ratio of 10% is used for the 
subsonic case. The relevant input parameters are shown in Figure 4.5.1. For the upper and 
lower boundaries of the Channel, the solid-wall boundary condition for invicid calculation, 

1. e. zero mass flex through the surface was used. The inflow boundary is on the left hand 
side. For subsonic inflow conditions, total pressure and total enthalpy are held constant. 
For the subsonic outflow conditions, we define only the back pressure (through the initial 
condition PIN) and extrapolate other flow parameters from the last grid points in the 
calculation domain. Figure 4.5.2 and 4.5.3 present isomach lines, as well as the Mach 
number distributions along the top and bottom walls of the steady invicid subsonic flows 
solution at inlet Ma=0.5. The symmetric Mach number contours about the center chord 
indicates the accuracy of the numerical solution. This 65x17 solution takes 9.67 CPU 
seconds for 50 time steps and uses 0.24 MW sore memory on a CRAY XMP/24. 

6. Transonic Row over a Circular Bump. 

This case uses the same geometry and grid points as the previous case except the 
inlet conditions. The relevant parameters are shown in Figure 4.6.1. For this transonic 
flow calculation with inlet Ma=0.675, a supersonic pocket appears in the solution and is 
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terminated by a shocks as seen from Figure 4.6.2 and 4.6.3 in which isomach contours 
and surface Mach profiles are presented. The shock is captured vividly with two transition 
points at the bottom wall. The shock is located at a distance 73% of the chord and the 
maximum value of Mach value is 1.34. These agree well with the values obtained by 
[23,24] for the same problem. The solutions reach steady state in 120 time steps with 
At=0.002 and take 23.14 CPU seconds. 

7. Superso nic Flow over a Circular Bump. 

This case involves a supersonic flow with inlet Ma=1.65 past a 4% bump using a 
99 x 33 grid. Supersonic inflow conditions requires specifications of all flow variables. For 
the supersonic outflow, all flow variables at the downstream boundary are extrapolated 
from the calculation domain points. These are done in the MAST code and are not directly 
reflected in the input file as shown in Figure 4.7.1. The calculated isomachs and surface 
Machs are shown in Figure 4.7.2 and 4.7.3 and compare well with the results of [25] 
using the second order Godunov method. A converged solution was obtained in 150 time 
steps with 51 CPU seconds and 0.46 MW core memory. 

8. Supersonic Flow Past a Blunt Body. 

To show the shock capturing capability involving real gases, two calculations were 
performed for hypersonic viscous flows, with free stream Mach number of 10.0 past a 
spherical nose of an axi-symmetric blunt body. The Reynolds number is 50,000 and the 
free stream temperature is assumed to be 300K. The standard k-e model with wall 

functions was used for the calculation. The input checks for the ideal gas case is shown in 
Figure 4.8.1 and the calculated isomachs are shown in Figure 4.8.2. Figure 4.8.3 and 
4.8.4 show the calculated temperature and Mach number along the stagnation line in front 
of the spherical nose. 

In addition to the ideal gas case, an equilibrium chemistry package based on Gibbs 
free energy minimization concept has been incorporated into the present method for real gas 
effects. An adiabatic wall boundary condition was employed for the solutions of energy 
equations for both cases. In the chemical equilibrium calculation, five species including O 2 , 
N 2 , NO, N and O are involved. Concentrations of these species are calculated based on 
minimization of Gibbs free energy. Input decks and calculated results are shown in Figures 
4.8.5, 4.8.6 and 4.8.7. This calculation was obtained by central different scheme with 
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50% upwind damping. To see the effect of damping, another calculation was run with 75% 
upwind damping (OMGPHI = 0.75). The results are shown in Figures 4.8.8 and 4.8.9. It 
can be seen that the shock is shaper for the 50% upwind damping case. However, the 
isomachs indicate oscillatory behavior for the lower, damping case. The CPU time 
requirements for this case are 101 seconds and 555 seconds for the ideal gas case and the 
equilibrium chemistry case respectively. The core memories are 0.46 MW and 0.58 MW 
respectively. 

9. Non-Evaporating Solid-Sprav 

The solid-cone spray measurement of Hiroyasu and Kadato [26] have used to 
validate the present dense spray model. Liquid fuel is injected through a single hole nozzle 
into constant pressure, room-temperature nitrogen. Spray tip penetration and drop sizes 
were measured from photographs of the backlighted spray. The nozzle diameter was 0.3 
mm and the present computations used tetradecane for the liquid fuel (the experiments used 
a diesel fuel oil with physical properties close to tetradecane). Atomization is introduced by 
injecting drops which have a characteristic size equal to the nozzle exit diameter (0.3 mm). 

The corresponding input data for P=l.l MPa are shown in Figure 4.9.1. A 
computational domain of 20 mm in radius and 120 mm in length was discretized by a 25 
radial and 45 axial grid. The mesh spacing was nonuniform with refinement on the 
centerline and close to the injector. The number of computational parcels was between 1000 
and 3000, and the number was varying with the back pressure. The present numerical 
results did not change appreciably when this parcel number was varied. The initial turbulent 
quantities were assumed as the small values (k = 1 x 10" 3 /n 2 /s 2 , e = 4 x 10~ 4 /ri 2 /s 3 ). The 
numerical results were insensitive to these initial values. 

Figure 4.9.2 and 4.9.3 shows the spray parcel distribution, and the predicted and 
measured spray tip penetration versus time. It can be seen that there is reasonably good 
agreement between the prediction and the measurement. A detailed discussions can be 
found in Refs.[27, 28]. 600 steps of this calculation take 480 CPU seconds and 0.42 MW 
core memory. 

10. Non-Evaporating Hollow-Cone Snrav 
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The hollow-cone spray tip penetration data of Shearer and Groff [29] have been 
used for the model validation. In the experiment, the liquid is injected into quiescent room- 
temperature nitrogen at P =550 k Pa. The corresponding input data are shown in Figure 
4.10.1. The numerical timestep is used as 2.5 us and about 2000 spray parcels are used in 
the computation. The experimental spray cone angle is 60 degrees, and the flow rate 
0.0165 mL/injection with four pulses, each of duration about 0.58 ms. The computational 
injection velocity are set to the experimental spray tip velocity ( 60 /d/s ) measured from the 
movie pictures in the early stage of the injection. 

Figure 4.10.2 and 4.10.3 shows the spray parcel distribution and the velocity 
vectors, and the predicted and measured spray tip penetration versus time. The velocity 
vectors show the presence of a vortex near the head of the spray, which curls the spray tip 
toward the outside of spray. A substantial region of strong inward flow in the center of the 
comer near the injector is also observed. The numerical results indicate that turbulence has a 
relatively small effect on penetration in a hollow-cone spray because the radial spreading of 
the spray dominates the spreading due to turbulence. It is seen that the predictions 
reasonably well agree with the experimental tip penetration. A detailed discussions can be 
found in Refs.[27, 28]. 528 steps of this calculation take 1020 CPU seconds and 0.67 MW 
core memory. 

11. Hollow-Cone Snrav Flames 

The present numerical model for the multi-phase turbulent reacting flows has been 
tested by applying it to predict the local flow properties in a axisymmetric, confined, 
swirling spray-combusting flows. Experimental data for temperature, axial and tangential 
velocity components were obtained from measurement of Khalil et. al. (30). The inlet 
conditions and the initial droplet size distribution used in the present computations are 
described in Ref.[31]. The corresponding input data are shown in Figure 4.11.1. Liquid 
kerosene was used as fuel and the air/fuel mass ratio was fixed at 20.17. 

In the present study, two swirling numbers (S = 0.72 and 1.98) were considered to 
investigate the influence of swirl on the droplet evaporation and burning characteristics. 
Figure 4.11.2, 4.11.3, and 4.11.4 show the general flow pattern such as the predicted 
droplet trajectories, velocity vectors, and temperature contours of two swirl cases. In the 
lower swirl case (S = 0.72), large portion of droplets survive in the central recirculation 
zone and continue to evaporate in the far downstream region. In the high swirl case (S = 
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1.98), most of small droplets are trapped in the recirculation zone and evaporate there, 
producing intensive burning and high temperature in this region. With increasing swirl, the 
droplet spreading increases due to the droplet dispersion and the increased particle 
centrifugal force term. In addition, the larger central recirculation zone corresponding to the 
higher inlet swirl is contributed to recirculate more hot combustion gas from downstream 
and to increase the temperature at near inlet regions. In this combusting dilute spray case, 
the present numerical procedure correctly predicts the general features of spray-combustion 
flows and yields the qualitative agreement with experimental data. The detailed discussions 
can be found from Refs.[27, 31]. The solutions reach the steady state in 1050 time steps 
with 2900 CPU seconds and 0.86 MW core memory. 
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CLOSURE 


In this report, we have briefly described the methodologies and some testings of the 
current version of MAST code. This code was written as an advanced research-oriented 
computer program, in contrast to a "black-box" production code. It is advisable that the 
user should have some knowledge and experience with computational fluid dynamics and 
pressure-based methodologies. It is necessary that the user become familiar with all parts 
of the program before attempting to use or modify it 

Although MAST has gone through a developing and debugging phase, by no 
means it is a fully mature program. More advanced physical submodels, as well as new 
ideas of numerics will be continuously coded and tested. Documentations of further 
developments and improvements of MAST will be published on a timely basis as soon as 
their validations and testings are completed. We would also liked to ask users to contact us 
if any anomalies were found while working with MAST. 

The research and development that went into the MAST code was supported by the 
NASA Mashall Space Flight Center under NASA Grants NAG8-092 and partially by 
NAG8 128. The authors are indebted to all of the people at NASA-MSFC that provided 
support, guidance, and assistance during the course of this study. 
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Table 2 . Source Terms S 0 and S* (continued) 
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TABLE 3-1 Descriptions of Subroutine 


MAIN 


AUXSUB (NCONTL) 


CADPXY (NCONTL) 
CALCAP 

CALCTB (NCONTL) 
CALCVU (NCONTL) 

CGSOL V (NCONTL), 
MATRIX, and DOTSUM 


Program MAIN handles the main control and sequence of 
prediction-corrector in each time step. It first calls routines 
INl l lL and TRANSF for setting up initial values and coefficients 
of grid transformation and then enters the main time-marching 
loop. Inside this loop, the predictor-corrector sequence is carried 
out. After the marching process is finished, MAIN also handles 
the writing of the solution to disks by calling OUTPUT. 

This is an auxiliary subroutine to cany out further some link 
coefficients for pressure increment equations (25) and (29). 
NCONTL = 1 and 2 are called from momentum equation to 
calculate phase values of H() in equation (24) and (27). NCONTL 
=3 and 4 are used for calculations of boundary pressure values after 
1st corrector and second corrector step. NCONTL =5 is used again 
for boundary pressures after all predictor- corrector steps are done. 

This is a utility subroutine to calculate pressure gradients for non- 
orthogonal grid. 

This subroutine computes link coefficients for all interior 
pressure-increment quantities except boundary points which are 
done by calling UPDTCF(2). 

This subroutiine computes corrector step of turbulent quantities k 
and e using explicit schemes. Boundary points are calculated 
through calling UPDTBC (4,4) and UPDTBC (4,5). 

This subroutine calculates both velocity components through 1st 
corrector step of (24) and 2nd corrector step of (27) using explicit 
expressions. NCONTL = 1 represents 1st corrector and NCONTL = 
2 represents 2nd corrector step. 

These subroutines activate the Conjugate gradient solver 
solution procedure for algebraic equations derived from the 
discretization procedure. NCONTL = 0 is used for pressure 
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COEFIS 


EXAMPL 


GEOMBC 


INITEL 


MODIFY 


NEWS TP 


SORCEM (NCONTL), 


increment finite difference equation solver and NCONTL = 1 is 
for the rest of the variables. 

This subroutine calculates link coefficients Ap and A n b for all 
transport equations which have the common form Apfp= SA n bfnb+ 
S u except two pressure incremental equations (25) and (29). 
Depending on the finite difference form used for discretization of 
convection terms, the link coefficients are constructed differently. It 
also calls UPDTCF for link coefficients of boundary points and 
MODCOF as well as MODSOR if modifications of these coefficients 
and source terms,which have the form Spfp , are required. 

Several tested example cases are supplied in this subroutine for the 
user to familiarize the current code. 

This subroutine calculates direction cosines and nonnal distance of 
the first interior points to the solid boundaries. 

This subroutine initializes most of the variables and arrays and reads 
in the input namelist to obtain the user-specified parameters. It calls 
KEYREAD to read in the key parameters (which will be described in 
KEYREAD). It contains some utilities such as reading in the restart 
file, grid systems if grids were generated by the user through other 
generation codes. Any blockage of calculation domain will also be 
done in this subroutine thus it calls GEOMBC to modify interior 
points. 

This subroutine gives user the freedom to create or modify boundary 
conditions, source terms and coefficients that are not included in the 
current version of the MAST code. 

The NEWSTP subroutine puts the most current calculated variables 
(i.e. variables at t n+ *) to the old time step before continuing the next 
time step cycle. 

These subroutines set up source terms for the 
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SORCEP (NCONTL), 
SORCSW, SORCTH, 
SORCET (NCONTL) 


SPRAY 


SYSTEM 


momentum eequation (NCONTL = 1 for u, 2 for v), 
pressre correction (NCONTL =1 for 1st 

corrector, 2 for second corrector), swirl velocity, thermal energy 
equation and turbulence quantities (NCONTL = 1 for k and 2 for £ ) 
respectively. 

This subroutine calculates spray droplets dynamics based on a 
lagrangian tracking scheme. 

Subroutine SYSTEM is the main control routine of the discretization 
method used in the MAST code. It composes of various ENTRY'S 
for setting up coefficients and sources of finite difference schemes, 
for solving the nonlinear algebraic equations and updating boundary 
conditions. For example, entry PREDIT solves equation (22) for 
velocity components in x direction (NF =1) and y direction (NF = 
2). It calls COEFIS to set up pressure terms, cross derivative 
diffusion terms and extra terms (such as swirl source terms and 
rotations if desired). 

The calling to CGSOLV is to utilize the conjugate gradient method 
for solveing the resulting discretized algebraic equations. 

Entry FSTCRT solves the incremental pressure (equation (25)) and 
incremental velocity equation (equation(24)). Note that equation (24) 
is explicit after pressure increment being solved thus by simply 
calling CALCVU is sufficient. Entry OTHCRT is deviced for 
second corrector step or other corrector steps if high order time 
accurate discretization in temporal domain were used. Entries 
THERML and SWIRL are used for thermal energy equation in heat 
transfer prolems and Swirling Velocity Component (w). Entries 
PREDTB and FCRTTB are predictor and corrector procedures for 
turbulence quantities such as kinetic energy equation (k) and its 
dissipation rate ( £ ) equation if the k- £ two equation model were 
used for turbulent flow calculation. It should be noted that the 
predictor step for turbrlent quantities is implicit and the corrector 
step is explicit 
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TRANSF 


This subroutine carries out grid transformation coefficients, and 
Jacobians by a central difference scheme. 


UPDTBC 

(NCONTL, NTUB) 


UPDTCF (NCONTL) 


This subroutine updates various boundary conditions 
for fluid flow and heat transfer problems. When NCONTL is set to 
be 1, it updates outflow boundary conditions based on mass 
conservation balance. 


NCONTL = 2 represents symmetric = 0 ) or slip boundary 

oX, 

conditions. NCONTL = 3 sets up cyclic boundary conditions. If 
wall function options were used in the turbulent flow calculation, 
NCONTL should be set to 4. Under this option, NTUB = 1 
recalculates an effective viscosity at the first point away from the 
wall. NTUB = 2 and 3 re-establish new source terms for k and £ 
equation in the first corrector step. NTUB = 4 and 5 create source 
terms for explicit k- e equation corrector step. NCONTL = 5 is 
used for adabatic boundary condition in thermal energy equation. 


This subroutine is used for updating the link coefficients for 
boundary points specifically for body-fitted coordinates when the 
grid lines do not intersect boundary orthogonally. NCONTL = 1 
activates this action for all variables except pressure. NCONTL = 2 
is used for pressure. 
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BEGIN 
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CONTROL NCRT 2 

INCOMP OMGD 0 OMGF 0.4 PHI -1.000 OMGPHI 0.00 


IMON 4 JMON 3 MONU 

ERRCG 

l.E-1 ERRM 1 

. E-4 




; RESTART 









GRID NX 

51 NY 51 







XDIR 

1ST 1 

IEND 

26 DST 0. 

0 

DEND 

0.5 

EXP 

1.50 

XDIR 

1ST 26 

IEND 

51 DST 0. 

5 

DEND 

1.0 

EXP 

-1.50 

YDIR 

1ST 1 

IEND 

26 DST 0. 

0 

DEND 

0.5 

EXP 

1.50 

YDIR 

1ST 26 

IEND 

51 DST 0. 

5 

DEND 

1.0 

EXP 

-1.50 

BOUND 









1ST 1 

I END 

51 JST 

51 JEND 

51 

WALL 

U 1 

♦ 


1ST 1 

IEND 

51 JST 

1 JEND 

1 

WALL 




1ST 51 

I END 

51 JST 

1 JEND 

51 

WALL 




1ST 1 

IEND 

1 JST 

1 JEND 

51 

WALL 




PROPERTY 

VISCOS 

1.0E-3 







SOLV 

U V 

P 







RUN 

DT 1.0 

NSTEP 150 NPR1 

. 10 




ENDJOB 










Figure 4.1.1 Input deck for the driven cavity flow 



X 


Figure 4.1.2 Grids for the driven cavity flow 
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CONVERGENCE HISTORY FOR DRIVEN CAVITY FLOW (Re=1000) 


STEP= 

STEP= 

STEP= 

STEP= 

STEP= 

STEP= 

STEP= 

STEP= 

STEP= 

STEP= 

STEP= 

STEP= 

STEP= 

STEP= 


10 TIME=1 . 000E+01 
20 TIME=2 . OOOE+Ol 
30 TIME=3. 000E+01 
40 TIME=4. 000E+01 
50 TIME=5. 000E+01 
60 TIME=6. OOOE+Ol 
70 TIME=7. OOOE+Ol 
80 TIME=8. OOOE+Ol 
90 TIME=9. OOOE+Ol 
100 TIME=1. 000E+02 
110 TIME=1. 100E+02 
120 TIME=1 . 200E+02 
130 TIME=1 • 3 00E+02 
140 TIME=1 . 400E+02 


FMON=-4 . 142E-05 NF= 
FMON= 6 . 257E-05 NF= 
FMON= 1. 794E-04 NF= 
FMON= 2 . 568E-04 NF= 
FMON= 2 . 766E-04 NF= 
FMON= 2 . 833E-04 NF= 
FMON= 2 . 878E-04 NF= 
FMON= 2.896E-04 NF= 
FMON= 2 . 888E-04 NF= 
FMON= 2.893E-04 NF= 
FMON= 2.895E-04 NF= 
FMON= 2 . 898E-04 NF= 
FMON= 2 . 898E-04 NF= 
FMON= 2.901E-04 NF= 


2 FERRM=1.023E-01 
2 FERRM=3 . 528E-02 
2 FERRM=1.546E-02 
2 FERRM=5 .97 1E-0 3 

1 FERRM=2 . 384E-03 

2 FERRM=1 . 424E-03 
2 FERRM=9 . 3 68E-04 
2 FERRM=6.548E-04 
2 FERRM=4 . 813E-04 
2 FERRM=3 . 512E-04 
2 FERRM=2 . 619E-04 
2 FERRM*=1 . 998E-04 
2 FERRM=1 . 490E-04 
2 FERRM=1 . 120E-04 


Figure 4.1.3 Sample monitors for the driven cavity flow 



1 ' 1 1 1 1 i 1 i 1 i 

0 02 0.4 06 OB 1 


X 

Figure 4.1.4 Computed streamlines for 
the driven cavity flow 
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Figure 4.2.1 Grid system for a 2-D flow 
over a circular cylinder 


CONTROL NCRT 2 

INCOMP OMGD 0. PHI -1.0 OMGPHI 0.00 OMGF 1.0 
IMON 4 JMON 3 MONU ERRCG l.E-2 ERRM l.E-8 
.•RESTART 

GRID NX 41 NY 41 

;GRID SET IN SUBROUTINE EXAMPLE 

BOUND 


1ST 

1 

I END 

41 

JST 

1 

JEND 

1 

WALL 

1ST 

1 

I END 

41 

JST 

41 

JEND 

41 

WALL U 1 . 

1ST 

1 

I END 

1 

JST 

1 

JEND 

41 

CYCLE 


PROPERTY VISCOS 5.0E-3 ;RENOLDS NUMBER=200 
SOLV U V P 

RUN DT 0.19 NSTEP 500 NPR1 10 NEX 2 

; DT 0.02625 NSTEP 50 NPR1 10 NEX 2 

ENDJOB 

Figure 4.2.2 Input deck for the cylinder flow 


39 





- SI- 




t= 100.25 

Figure 4.2.3 Streamline Patterns at different times within 
a shadding cycle 
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; BACKWARD FACING STEP, LAMINAR FLOW — 

CONTROL 

INCOMP OMGD 0 PHI -1.00 OMGPHI 0.0 ERRCG l.E-1 ERRM l.E-4 
IMON 17 JMON 4 MONU OMGF 0.2 
NCRT 3 
/RESTART 


GRID 

NX 

81 NY 

51 








XDIR 

1ST 

1 

IEND 

81 

DST 0 

. 0 

DEND 24. 

EXP 

+ 1.00 

YDIR 

1ST 

1 

IEND 

51 

DST 0 

. 0 

DEND 1 . 0 

EXP 

+ 1.00 

BOUND 











1ST 

1 

IEND 

11 

JST 

1 

JEND 

26 

BLOCK 



1ST 

1 

IEND 

1 

JST 

27 

JEND 

51 

INLET U 

1 


1ST 

81 

IEND 

81 

JST 

1 

JEND 

51 

OUTLET 



1ST 

12 

IEND 

81 

JST 

1 

JEND 

1 

WALL 



1ST 

1 

IEND 

81 

JST 

51 

JEND 

51 

WALL 




PROPERTY VISCOS 1.25E-3 DENGAS 1. ; RENOLDS NUMBER=800 

SOLV U V P 

RUN DT 1.25 NSTEP 350 NPR1 10 NEX 1 

ENDJOB 

Figure 4.3.1 Input deck for the Laminar Backward Step Flow 



X 

Figure 4.3.2 Contour of Streamlines for the Laminar Backward 
fhcing Step Flow 
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CONTROL NCRT 3 

INCOMP OMGD 0 PHI -1.0 OMGPHI .30 ERRCG l.E-1 ERRM l.E-4 
; RESTART 

IMON 15 JMON 4 MONU 

GRID NX 59 NY 35 ; XLEN 18. YLEN 1. 


XDIR 

1ST 

1 

IEND 

11 

DST 

0 . 

DEND 

1.5 

EXP 

-1.5 

XDIR 

1ST 

11 

I END 

40 

DST 

1.5 

DEND 

4.5 

EXP 

1.2 

XDIR 

1ST 

40 

IEND 

59 

DST 

4.5 

DEND 

18. 

EXP 

1.5 

YDIR 

1ST 

1 

IEND 

15 

DST 

0 . 

DEND 

0.3333 

EXP 

1.0 

YDIR 

1ST 

15 

IEND 

25 

DST 

0.3333 

DEND 

0.6666 

EXP 

1.5 

YDIR 

1ST 

25 

IEND 

35 

DST 

0.6666 

DEND 

1.0 

EXP 

-1.5 


BOUND 

1ST 

1 

IEND 

59 

JST 

1 

JEND 

1 

WALL 

1ST 

1 

IEND 

59 

JST 

35 

JEND 

35 

WALL 

1ST 

1 

IEND 

1 

JST 

16 

JEND 

35 

INLET U 1 

1ST 

59 

IEND 

59 

JST 

1 

JEND 

35 

OUTLET 

1ST 

1 

IEND 

11 

JST 

1 

JEND 

15 

BLOCK 

PROPERTY VISCOS 1. 

000E-5 

; RENOLDS NUMBER=1.E5 


TURBULEN TKIN 0.003 SCALE 0.33 

SOLV U V P TK TE ; TURBULENCE MODEL ACTIVATED 

RUN DT 0.6 NSTEP 90 NPR1 10 

ENDJOB 

Figure 4.4.1 Input deck for the turbulent backward- facing step flow 
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Figure 4.4.2 Grid system for the backward- facing step flow 







17 


JST 


CONTROL 

COMPRES OMGD 0.0 PHI +0.33 OMGPHI 0.00 
NCRT 3 ERRCG l.E-1 ERRM 1 . E-4 IMON 4 
; RESTART 

GRID NX 65 NY 
BOUND 

1ST 1 IEND 
U 169.37754 
1ST 65 IEND 
1ST 1 IEND 
1ST 1 IEND 
PROPERTY VIS COS 
UIN 169.37754 
SOLV U V 


TEMP 285.7143 


1 

1 

17 


65 JST 
65 JST 
65 JST 
O.E-3 

TIN 285.7143 
P TEMP 


JEND 

P 

JEND 

JEND 

JEND 


OMGT 0.5 OMGF 0.1 
JMON 3 MONU 


17 INLET 
1.0E5 ;Min=0 . 5 

17 OUTLET 
1 SLIP Q 0. 

17 SYMMETRY Q 0. 


PSTAG 1 . 18621E5 TSTAG 300 PIN ] 


RUN 

END JOB 


DT 3.0E-3 


NSTEP 50 NPR1 10 NPR2 600 


NEX 8 


Figure 4.5.1 Input deck for the subsonic bump flows 



Figure 4.5.2 IsoMach line contours of the subsonic bump flow 


o 



Figure 4.5.3 Mach number distributions along the walls 


. 0E5 
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CONTROL - 

COMPRES OMGD 1.0 PHI +0.90 OMGPHI 0.15 OMGT 0.8 OMGF 0.10 
NCRT 3 ERRCG l.E-1 ERRM 1 . E-4 IMON 4 JMON 3 MONU 
; RESTART 


GRID 

BOUND 

NX 

65 NY 

17 





1ST 

1 

IEND 

1 

JST 1 

JEND 

17 

INLET 

U 221.20149 

TEMP 274.9456 P 

1 . 

0E5 ;Min=0. 675 

1ST 

65 

IEND 

65 

JST 1 

JEND 

17 

OUTLET 

1ST 

1 

IEND 

65 

JST 1 

JEND 

1 

SLIP Q 0. 

1ST 

1 

IEND 

65 

JST 17 

JEND 

17 

SYMMETRY Q 0. 

PROPERTY 

VISCOS 

0.E- 

-3 




UIN 

221 

.20149 

TIN 

274.9456 

PS TAG 

1 . 

35694E5 TSTAG 300 


SOLV U V P TEMP 

RUN DT 2.0E-3 NSTEP 120 NPR1 10 NPR2 600 NEX 8 

ENDJOB 

Figure 4.6.1 Input Deck for the Transonic Bump Flow 



Figure 4.6.2 IsoMach Contours for the Transonic Bump Flow 


iO 



Figure 4.6.3 Surface Mach Numbers Along the Walls 


. 0E5 
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CONTROL 

COMPRES OMGD 1.0 PHI +0.33 OMGPHI 0.10 OMGT 0.5 OMGF 0.5 
NCRT 3 ERRCG l.E-1 ERRM l.E-4 IMON 4 JMON 3 MONU 
; RESTART 

GRID NX 99 NY 33 
BOUND 

1ST 1 I END 1 JST 1 JEND 33 INLET 

U 460.86153 TEMP 194.2376 P 1.0E5 ;Min=1.65 DT 0.3E-3 OMGT 0.5 

1ST 99 I END 99 JST 1 JEND 33 OUTLET 

1ST 1 IEND 99 JST 1 JEND 1 SLIP Q 0. 

1ST 1 IEND 99 JST 33 JEND 33 SYMMETRY Q 0. 

PROPERTY VISCOS O.E-3 

UIN 460.86153 TIN 194.2376 PSTAG 4.57886E5 TSTAG 300 PIN 1.0E5 
SOLV U V P TEMP 

RUN DT 0.3E-3 NSTEP 150 NPR1 10 NPR2 600 NEX 9 

ENDJOB 

Figure 4.7.1 Input Deck for the Supersonic Bump Flow 



-15 -1 -05 6 05 i 15 


X 

Figure 4.7.2 IsoMach Contours of the Supersonic Pump Flow 



-is -1 -05 0 05 1 15 

X 


Figure 4.7.3 Surface Mach Numbers Along the Walls 
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CONTROL NCRT 3 

COMPRES OMGD 1 PHI +1.0 OMGPHI 0.50 OMGT 0.3 OMGF .00 
;FOR EQUILIBRIUM 0MGPHI=0.75 
IMON 4 JMON 3 MONU ERRCG l.E-1 ERRM l.E-4 
? RESTART 

GRID NX 81 NY 41 AXISYM 
BOUND 


1ST 

1 

I END 

1 JST 

1 

JEND 

41 

SYMMETRY 

Q 0 

1ST 

81 

IEND 

81 JST 

1 

JEND 

41 

OUTLET 


1ST 

1 

I END 

81 JST 

1 

JEND 

1 

WALL U 0. 

V 0. 

1ST 

1 

IEND 

81 JST 

1 

JEND 

1 

SLIP Q 

0 

1ST 

1 

IEND 

81 JST 

41 

JEND 

41 

INLET 



U 3471.269240 TEMP 300. P 1.658E4 TK 10000. TE 30000000. 

PROPERTY VISCOS 2.00E-5 

UIN 3471.269240 TIN 300. PSTAG 7.1510109E8 TSTAG 6300. PIN 1.658E4 
TURBULEN 

TKIN 10000. TEIN 30000000. ; SCALE 1.0 
SOLV U V P TEMP TK TE 

RUN DT 1.5E-4 NSTEP 50 NPR1 10 NEX 13 

ENDJOB 


Figure 4.8.1 Input Deck for the Hypersonic Flow Past a 
Blunt Body- Ideal Gas 



X 

Figure 4.8.2 Mach Number Contours for the Blumt Body Flow 
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MACH 


-4 -35 -3 

X 

Figure 4.8.3 Temperature Profile Along the Stagnation Line 



Figure 4.8.4 Mach Number Along the Stagnation Line 
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CONTROL NCRT 3 

COMPRES OMGD 1 PHI +1.0 OMGPHI 0.50 OMGT 0.3 OMGF .00 
; FOR EQUILIBRIUM OMGPHI=0.75 
IMON 4 JMON 3 MONU ERRCG l.E-1 ERRM l.E-4 
; RESTART 

GRID NX 81 NY 41 AXISYM 
BOUND 


1ST 

1 

IEND 

1 JST 

1 

JEND 

41 

SYMMETRY Q 0 

1ST 

81 

IEND 

81 JST 

1 

JEND 

41 

OUTLET 


1ST 

1 

IEND 

81 JST 

1 

JEND 

1 

WALL U 

• 

o 

> 

• 

o 

1ST 

1 

IEND 

81 JST 

1 

JEND 

1 

SLIP 

Q 0 

1ST 

1 

IEND 

81 JST 

41 

JEND 

41 

INLET 


U 3471. 

269240 

TEMP 300. 

P 

1.658E4 

TK 

10000. 

TE 30000000 


PROPERTY VISCOS 2.00E-5 

UIN 3471.269240 TIN 300. PSTAG 7.1510109E8 TSTAG 6300. PIN 1.658E4 
TURBULEN 

TKIN 10000. TEIN 30000000. ; SCALE 1.0 

REACTION 

NSPE 5 NELE 2 NFROS 2 HP ; FROSEN TP HP 
ELEMCONS 1 3.76 

CHEMFORM ; INPUT DATA FORMAT ( 4X , A8 , 8F2 . 0 ) 

02 2 0 

N2 0 2 

NO 11 

N 0 1 

0 10 

SOLV CONHT EQLM U V P TEMP TK TE 
RUN DT 1.5E-4 NSTEP 150 NPR1 10 NEX 13 

ENDJOB 

Figure 4.8.5 Input deck for hypersonic flow-equilibrium air 





Figure 4.8.6 Contour of Mach number 
50% damping 
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Figure 4.8.7 Stagnation Line Temperature for Ideal Gas Air 
and Equilibrium air, 501 Damping 



X 

Figure 4.8.8 Mach Number Contours, 751 Damping 
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75% damping 
50% damping 



-4 -35 -3 

X 

Figure 4.8.9 Effects of Upwing Damping on Stagnation 
Line Temperature Profiles 
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; SOLID-CONE SPRAY 

CONTROL COMPRES 

OMGD 1 PHI -1.00 OMGPHI 1.000 ERRCG l.E-3 ERRM +1.E-7 
NCRT 2 OMGF 1.0 IMON 10 JMON 10 MONU 
/RESTART 

GRID NX 45 NY 25 AXISYM 


XDIR 

1ST 

1 

IEND 

2 

DST 

-l.E-3 

DEND 

0.000 

EXP 

1.0 

XDIR 

1ST 

2 

IEND 

45 

DST 

0.00 

DEND 

0.12 

EXP 

1.3 

YDIR 

1ST 

1 

IEND 

10 

DST 

0 . 

DEND 

5.E-3 

EXP 

1.0 

YDIR 

1ST 

10 

IEND 

25 

DST 

5.E-3 

DEND 

0.020 

EXP 

1.3 


BOUND 


1ST 

1 

IEND 

45 

JST 

1 

JEND 

1 

SYMMETRY 


1ST 

1 

IEND 

45 

JST 

25 

JEND 

25 

OUTLET P 

1.1E+6 

1ST 

1 

IEND 

1 

JST 

1 

JEND 

25 

WALL 


1ST 

45 

IEND 

45 

JST 

1 

JEND 

25 

OUTLET P 

1.1E+6 


PROPERTY VISCOS 1.7E-5 DENGAS 12.36 UIN 0.00 TIN 298 PIN 1.1E+6 
YNIN 1.0 YFIN 0.0 YOIN 0.0 
TURBULEN TKIN 1.0E-3 TEIN 4 . OE-4 

SPRAY SMR 150. E-6 CONST DENPT 840 TKESW BREAKUP COLIDE 

CONE2 10 CONE1 0 MODUL 0 
1ST 2 IEND 2 JST 2 JEND 2 NPTS 4 
VINJ 115.8 FLOWP 6.876E-3 
END 

EVAP TEMP 298 TABEL 

SOLV U V P TK TE PATC 

RUN DT 10.0E-6 NSTEP 250 NPR1 10 NPR2 50 

ENDJOB 


Figure 4.9.1 Input deck for solid-cone spray 
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; SOLID-CONE SPRAY 

CONTROL COMPRES 

OMGD 1 PHI -1.00 OMGPHI 1.000 ERRCG l.E-3 ERRM +1.E-7 
NCRT 2 OMGF 1.0 IMON 10 JMON 10 MONU 
; RESTART 

GRID NX 45 NY 25 AXISYM 


XDIR 

1ST 

1 

IEND 

2 

DST 

-l.E-3 

DEND 

0.000 

EXP 

1.0 

XDIR 

1ST 

2 

IEND 

45 

DST 

0.00 

DEND 

0.12 

EXP 

1.3 

YDIR 

1ST 

1 

IEND 

10 

DST 

0 . 

DEND 

5.E-3 

EXP 

1.0 

YDIR 

1ST 

10 

IEND 

25 

DST 

5.E-3 

DEND 

0.020 

EXP 

1.3 


BOUND 


1ST 

1 

IEND 

45 

JST 

1 

JEND 

1 

SYMMETRY 


1ST 

1 

IEND 

45 

JST 

25 

JEND 

25 

OUTLET P 

3 . OE+6 

1ST 

1 

IEND 

1 

JST 

1 

JEND 

25 

WALL 


1ST 

45 

IEND 

45 

JST 

1 

JEND 

25 

OUTLET P 

3 . OE+6 


PROPERTY VISCOS 1.7E-5 DENGAS 33.70 UIN 0.00 TIN 298 PIN 3 . OE+6 
YNIN 1.0 YFIN 0.0 YOIN 0.0 
TURBULEN TKIN 1.0E-3 TEIN 4 . 0E-4 

SPRAY SMR 150. E-6 CONST DENPT 840 TKESW BREAKUP COLIDE 

CONE 2 10 CONE1 0 MODUL 0 

1ST 2 IEND 2 JST 2 JEND 2 NPTS 5 
VINJ 102.5 FLOWP 6.090E-3 
END 

EVAP TEMP 298 TABEL 

SOLV U V P TK TE PATC 

RUN DT 18.0E-6 NSTEP 250 NPR1 10 NPR2 10 

ENDJOB 


; SOLID-CONE SPRAY 

CONTROL COMPRES 

OMGD 1 PHI -1.00 OMGPHI 1.000 ERRCG l.E-3 ERRM +1.E-7 
NCRT 2 OMGF 1.0 IMON 10 JMON 10 MONU 
; RESTART 

GRID NX 45 NY 25 AXISYM 

XDIR 1ST 1 IEND 2 DST -l.E-3 DEND 0.000 EXP 1.0 
XDIR 1ST 2 IEND 4 5 DST 0.00 DEND 0.12 EXP 1.3 

YDIR 1ST 1 IEND 10 DST 0. DEND 5.E-3 EXP 1.0 

YDIR 1ST 10 IEND 25 DST 5.E-3 DEND 0.020 EXP 1.3 

BOUND 


1ST 

1 

IEND 

45 

JST 

1 

JEND 

1 

SYMMETRY 


1ST 

1 

IEND 

45 

JST 

25 

JEND 

25 

OUTLET P 

5 . OE+6 

1ST 

1 

IEND 

1 

JST 

1 

JEND 

25 

WALL 


1ST 

45 

IEND 

45 

JST 

1 

JEND 

25 

OUTLET P 

5. OE+6 


PROPERTY VISCOS 1.7E-5 DENGAS 56.17 UIN 0.00 TIN 298 PIN 5. OE+6 
YNIN 1.0 YFIN 0.0 YOIN 0.0 
TURBULEN TKIN 1.0E-3 TEIN 4.0E-4 

SPRAY SMR 150. E-6 CONST DENPT 840 TKESW BREAKUP COLIDE 

CONE 2 10 CONE1 0 MODUL 0 

1ST 2 IEND 2 JST 2 JEND 2 NPTS 5 
VINJ 86.41 FLOWP 5.130E-3 
END 

EVAP TEMP 298 TABEL 

SOLV U V P TK TE PATC 

RUN DT 18.0E-6 NSTEP 250 NPR1 10 NPR2 10 

ENDJOB 
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Figure 4.9.2 Spray parcel distribution in a solid-cone spray (t=1.0ms) 
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; HOLLOW-CONE SPRAY 

CONTROL COMPRES 

OMGD 1 PHI -1.00 OMGPHI 1.000 ERRCG l.E-3 ERRM +1.E-7 


NCRT 2 OMGF 1. 

0 

IMON 

10 

JMON 

10 

MONU 




; RE START 










GRID NX 81 NY 

41 

AXISYM 







XDIR 1ST 

1 

IEND 

81 

DST 

0.00 

DEND 

0.040 

EXP 

1.0 

YDIR 1ST 

1 

IEND 

41 

DST 

0.00 

DEND 

0.020 

EXP 

1.0 


BOUND 


1ST 

1 

IEND 

5 

JST 

1 

JEND 

3 

BLOCK 



1ST 

6 

IEND 

81 

JST 

1 

JEND 

1 

SYMMETRY 


1ST 

1 

IEND 

1 

JST 

1 

JEND 

41 

WALL 



1ST 

1 

IEND 

81 

JST 

41 

JEND 

41 

OUTLET 

P 

0. 55E+6 

1ST 

81 

IEND 

81 

JST 

1 

JEND 

41 

OUTLET 

P 

0.55E+6 


PROPERTY VISCOS 1.7E-5 DENGAS 6.360 UIN 0.00 TIN 300 PIN 0.55E+6 
YNIN 1.0 YFIN 0.0 YOIN 0.0 
TURBULEN TKIN 1.0E-3 TEIN 4 . OE-4 

SPRAY SMR 30.E-6 CONST DENPT 840 TKESW BREAKUP COLIDE 

CONE1 30 CONE 2 30 MODUL 0 
1ST 6 IEND 6 JST 3 JEND 3 NPTS 5 
VINJ 60 FLOWP 5.980E-3 
END 

EVAP TEMP 300 TABEL 

SOLV U V P TK TE PATC 

RUN DT 5.E-6 NSTEP 200 NPR1 10 NPR2 50 

ENDJOB 


Figure 4.10.1 Input deck for hollow-cone spray 
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PENETRATION (mm) 



Figure 4.10.3 Spray tip penetration versus time 
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HOLLOW-CONE SPRAY FLAME 

HYBRID SWIRL IMON 4 JMON 3 MONU ERRCC 1 . E-l ERRM IE-8 
IPF 2 JPF 2 PREF 1.0568E+5 

RESTART 

GRID NX 61 NY 41 AXISYM 

XDIR 1ST 1 I END 61 DST 0.00 DEND 2.000 EXP 1.3 
YDIR 1ST 1 I END 41 DST 0 DEND 0.100 EXP 1.0 

BOUND 

1ST 1 I END 61 JST 1 JEND 1 SYMMETRY 

1ST 1 I END 61 JST 41 JEND 41 WALL Q 0 

1ST 1 I END 1 JST 1 JEND 8 WALL Q 0 

1ST 1 I END 1 JST 9 JEND 19 INLET U 15.0 OMEGA 22.3 

TK 1.125 YFU 0.0 YOX 0.233 G 0.0 
1ST 1 I END 1 JST 20 JEND 41 WALL Q 0 

1ST 61 I END 61 JST 1 JEND 41 OUTLET 

PROPERTY VISCOS 2.0E-5 TEMP 310 KGAS 0.026E+3 CPGAS 1000 
DENGAS 1.183 

TURBULEN CT1 1.44 CT2 1.92 CMU 0.09 SME 1.3 SMK 1.0 
TKIN 1.125 

SPRAY SMR +6.35E-5 DENPT 770 TEMP 310 EVAP 
MODUL 0 

1ST 1 I END 1 JST 1 JEND 1 NPTS 15 VINJ 11. FLOWP 4.89E-3 
SOLV U V P TK TE SW FU PATC G 

RUN DT 2.0E-4 NTIME 250 NPR1 10 NPR2 10 

ENDJOB 

Figure 4.11.1 Input deck for hollow-cone spray flame 
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Figure 4.11.2 Droplet trajectories in kerosene spray flame fields 
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Figure 4.11.3 Velocity vectors in kerosene spray flame fields 



APPENDIX A 


List of FORTRAN Symbols 


AE() 

AIRDIF 

AIRLA1 

AIRLA2 

AIRMU1 

AIRMU2 

AK() 

AKG 

AKPAS( ) 

ALF() 

AM() 

AMMR 

AMW( ) 

AN() 

ANE( ) 
ANW( ) 
AP() 

APO( ) 
AS() 

ASE( ) 
ASW( ) 
AW() 

BOO 

BIG 

BREAKUP 

CMU 

CONE1 

CONE2 

CP() 

CPG 
CPPAS( ) 
CSUBK 
CSUBMU 


discretization coefficient of east point 
molecular diffusivity of fuel vapor in ambient medium 
is AIRDIF*TEMP**EXPDIF / RHO 
molecular heat conduction coefficient of ambient medium 
is AIRLA1 *TEMP** 1.5/ (TEMP+AIRLA2) 
molecular viscosity of ambient medium 

is AIRMU 1 *TEMP** 1.5/ (TEMP+AIRMU2) 
heat conductivity k 
constant heat conductivity 
heat conductivity at patch 

declination angle in radian between grid line and boundary 
Mach number 
droplet mass median radius 
molecular weighting factor 
discretization coefficient of north point 
discretization coefficient of north east point 
discretization coefficient of north west point 
discretization coefficient of main point 
coefficient of CG solver 
discretization coefficient of south point 
discretization coefficient of south east point 
discretization coefficient of south west point 
discretization coefficient of west point 
coefficient of CG solver 
a very large real number 
switch of BREAKUP model(l - on, 0 - off) 
turbulence model constant 
fuel injector spray cone angle 
( CONE2>CONE 1 ) 
specific heat 
constant specific heat 
specific heat at patch 
constant of BREAKUP model 
constant of BREAKUP model 


CSX() 
CSY() 
CTX( ) 

CTY( ) 
CT1 
CT2 
D0() 
DELT( ) 
DEN( ) 
DENO( ) 


^/3x 

d^/ay 

dT|/dx 

dr\/dy 

turbulence model constant 
turbulence model constant 
coefficient of CG solver 

first grid cell size of non-uniform grid section 
density 

density of last time step 
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DENA( ) 
DENGAS 
DENPAS( ) 
DET( ) 
DLEND( ) 
DLH( ) 
DLST( ) 
DMAX 
DMTOT( ) 
DPX( ) 
DPY( ) 
DSIEP( ) 
DSX( ) 
DSY( ) 

DT 

DTEPS( ) 
DTKEP( ) 
DTMAX 
DTMIN 
DU() 

DVO 
EK( ) 

ELIQ( ) 

ENTHO( ) 

EPS() 

ERRCG 

ERRM 

EVAPM 

EVAPP 

EXPDIF 

EXPQ 

FERR( ) 

FEXP( ) 

FINP( ) 

FLOWP() 

FN( ) 

FP( ) 

GAM() 

GAMF 

GAMI 

GAMMA( ) 

GENT( ) 

GI() 

GJ() 

GMW 

GPI() 

GPJ() 

GX 

GY 


density of the time step before the last 
constant density 
density at patch 

Jacobian of coordinate tranformation 

end position of non-uniform grid section 

distance between boundary and first grid point from wall 

start position of non-uniform grid section 

control criterion in CG solver 

particle contribution to cell mass 

pressure gradient in X direction 

pressure gradient in Y direction 

particle contribution to cell specific internal energy 

distance between grid point in X direction 

distance between grid point in Y direction 

time step 

particle contribution to cell turbulence £ 
particle contribution to cell turbulence k 
maximum time step 
minimum time step 

connection coefficients of pressure correction and 
velocity correction 
a conversion of HK array 
a conversion of HLATO array 
initial cell enthalpy in EVAP 
e of turbulence model 

convergence criterion for conjugate gradient matrix solver 
convergence criterion for steady state solution 
particle evaporation mass in domain 
switch of evaporation model 
( see AIRD1F ) 

constant for particle size distribution 

relative error of each time step 

stretching factor of non-uniform grid 

interpolation function for boundary point 

injected particle mass flow rate 

array to store all the main variables 

temporary array 

T of differential equation 

constant for particle size distribution 

constant ratio of specific heat 

ratio of specific heat 

turbulence generation term 

pressure damping term used in scheme 

pressure damping term used in scheme 

molecular weighting factor 

pressure damping term used in face velocity 

pressure damping term used in face velocity 

gravity in X direction ( used in particle phase ) 

gravity in Y direction ( used in particle phase ) 
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GZ 

Hon 

HB( ) 

HK() 
HLATO( ) 
HPAS() 
HPTC( ) 
HPU( ) 
HPV( ) 
HSTAG 
HTFORM( ) 
IBC() 
IDBC( ) 
IDFI( ) 
IDSP( ) 
IEND( ) 

I LEND ( ) 
ILST( ) 
IMON 
INCOMP 
INEB( ) 
INX() 

INY() 
IPEND( ) 
IPST( ) 
IPTC( ) 

IST( ) 
JBC() 
JEND( ) 
JMON 
JPEND( ) 
JPST( ) 
JPTC( ) 
JST() 
KOLIDE 
KSIZE 

LDBC( ) 

LF 

LH 

LO 

LP 

LPATC 

LPFAC 

LPGEO 

LPT1 

LPT2 

LREST 

LS() 


gravity in Z direction ( used in particle phase ) 
coefficient of CG solver 
coefficient of CG solver 
enthalpy of species at T=100(N-1) 
latent heat of liquid 
enthalpy at patch 
for generating particle size 
operator for finite-difference representation of 
momentum equation 
total enthalpy 

heat of formation of species 

I index of boundary point 

type of boundary condition 

flag of variables specified at boundary 

identify species 

patch ending point 

end point of non-uniform grid section 

start point of non-uniform grid section 

monitor point I index 

switch of compressible ( 1 - incomp. 0 - comp. ) 

specify first grid point from wall 

relative I index of first grid point from wall 

relative J index of first grid point from wall 

particle injection position, ending 

particle injection position, starting 

particle index at cell 

patch start point in I direction 

J index of boundary point 

patch ending point in J direction 

J index of monitor point 

particle injection position.ending 

particle injection position, starting 

particle index at cell 

patch start point in J direction 

switch for particle collision model 

specify particle size distribution (0 - constant, 

1 - modified Rosin-Rammler, 2 - Chi-Squared) 
boundary condition type of patch 
fuel variable index in FN( ) 
index for enthalpy in FN( ) 
oxygen index in FN( ) 
index for pressure in FN( ) 
switch for particle tracking 
switch for face velocity printing 
switch for geometry parameter printing 
printing flag 
printing flag 
switch of restart 

index of variable need to be solved 
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LSTOP 

LTE 

LTK 

LU 

LV 

LW 

MBK 

MBL( ) 

MNP 

MODUL 

MPM 

MSP 

MSPI 

MVAP 

MX 

MY 

NAS 

NBC 

NCHEM 

NCGO 

NCGM 

NCRT 

NDIR( ) 

NEX 

NF 

NFL 

NFMAX 

NFMN 

NGRD 

NP 

NPC( ) 

NPR1 

NPR2 

NPTM 

NPTS( ) 

NSP 

NSTEP 

NSW 

NT 

NUMNOZ 

NX 

NX1 

NX2 

NY 

NY1 

NY2 

OMG() 

OMGD 

OMGF 


identification of stop of time marching 

index of turbulence model £ in FN( ) 

index of turbulence model k in FN( ) 

index of U velocity in FN( ) 

index of V velocity in FN( ) 

index of swirl velocity in FN( ) 

maximum dimension of patch related array 

identification of block point 

maximum dimension of particle arrays 

particle turbulence modulation models 

maximum dimension of boundary related arrays 

maximum number of species 

maximum dimension of 2-D species arrays 

maximum dimension of fuel property variables 

maximum dimension of I direction 

maximum dimension of J direction 

switch of axis symmetry 

number of boundary point 

switch of chemical reaction 

number of CG iteration used 

maximum CG solver iteration 

number of pressure correction steps 

direction of non-uniform grid 

index of example case 

index of variable in FN( ) 

switch of laminar flow 

maximum index of FN( ) 

index of monitor variable 

number of grid section 

number of patch 

particle number in a certain grid point 
output frequency 
output frequency 

current particle number in domain 

particle number injected per time step 

number of species to be solved 

maximum time marching step 

switch of swirl flow 

number of time marching steps 

number of nozzles 

number of point in I direction 

NX-1 

NX-2 

number of point in J direction 

NY-1 

NY-2 

rotation per second 

weighting factor for continuity density convection term 
face velocity interpolation weighting factor 
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OMGPHI 

OMGT 
OSCILO( ) 

P() 

POO 

pin 

PARTN( ) 
PCRIT 
PFO( ) 

PHI 

PI 

PI403R 

PIN 

PPAS( ) 

PRG 

PRT 

PSFA 

PSFB 

PSTAG 

PTEMP 

PVAP( ) 

QPASO 

R0() 

RIO 

RADP( ) 
RADPP( ) 
RCON( ) 
RCONI 
RELVEL( ) 
RERF( ) 
RHOO 
RHOI( ) 
RHOP 
RMW( ) 
RPMAX 
RU() 

RV() 

RW( ) 

sen 

SCALE 

SIE() 

SIXTH 

SMALL 

SMAX 

SME 

SMK 

SMR 


weighting factor of upwinding in Chakravarthy and 
Osher scheme 

supersonic temperature field relaxation 
initial droplet oscillation frequency at injector in 
BREAKUP 
pressure 

coefficient of CG solver 
coefficient of CG solver 

particle number in the computational parcel 

liquid phase critical pressure 

pressure of last time step 

parameter of Chakravarthy and Osher scheme 

n 

4/3*PI*RHOP 

initial pressure 

pressure at patch 

laminar Prandtl number 

turbulent Prandtl number 

constant for liquid phase vapor pressure 

constant for liquid phase vapor pressure 

stagnation pressure 

initial particle temperature 

liquid phase vapor pressure 

heat flux at patch 

coefficient of CG solver 

coefficient of CG solver 

particle size before evaporation 

particle size after evaporation 

gas constant 

initial gas constant 

particle-gas relative velocity 

inverse error function 

continuous phase density = DEN( ) 

initial species density 

particle density 

1/AMW 

maximum particle size 
particle contribution to cell momentum 
particle contribution to cell momentum 
particle contribution to cell momentum 
source term of discretization equation 
length scale for initial turbulence e 
enthalpy = FN(I,J,LH) 

1/6 

a very small real number 
maximum mass residual in continuity 
turbulence model constant 
turbulence model constant 
Sauter mean radius 
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SMT 

SP() 

SPDRAG( ) 

SPKT( ) 

SPMTIL( ) 

SSUM 

ST350 

STB 

STM 

SUVWO 

SWC() 

SWP( ) 

TOO 

TABEL 

TBN 

TCRIT 

TEIN 

TEMP( ) 

TEPAS( ) 

THIRD 

TIMAX 

TIME 

TIN 

TKE( ) 

TKESW 

TKIN 

TKPAS( ) 

TN( ) 

TOTCM( ) 
TOTH( ) 
TP() 
TPAS( ) 
TS() 
TSTAG 
TURBT( ) 
TWOTHD 
U() 

UCFO 
UCFO( ) 
UIN 
UP() 
UPAS( ) 
US2( ) 
UTRB( ) 
V() 

VAPM( ) 
VCF( ) 
VCFO( ) 
YIN 


turbulence model constant 
discretization coefficient of source term 
particle drag 

turbulence source term splitting 

explicit portion of contribution to mass diffusion 

summation of pressure correction equation source term 

calculation of particle surface tension 

calculation of particle surface tension 

calculation of particle surface tension 

particle coupling term for momentum exchange 

swirl velocity source term 

swirl velocity source term 

temperature of last time step 

flag for property using table or correlation 

liquid phase normal boiling temperature 

liquid phase critical temperature 

initial turbulence £ 

temperature = TN( ) 

turbulence e at patch 

1/3 

maximum time step 
current time marching 
initial temperature 
turbulence k 

flag for particle turbulent dispersion 
initial turbulence k 
turbulence k at patch 
continuous phase temperature 
total cell mass in EVAP 
total cell enthalpy in EVAP 
particle temperature 
temperature at patch 
intermediate step temperature 
stagnation temperature 

turbulence correlation time for particle diffusion 
2/3 

U velocity 

face velocity U 

last time step face velocity U 

initial U velocity 

particle velocity U 

U velocity at patch 

temporary array 

gas phase turbulent velocity fluctuation 

V velocity 

vapor mass in cell 

face velocity V 

last time step face velocity V 

initial V velocity 
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VINJ( ) 

particle injection velocity 

VISCOS 

viscosity 

VISLIQ( ) 

liquid phase viscosity 

VISM( ) 

laminar viscosity 

VIST( ) 

turbulent viscosity 

VOL( ) 

volume of cell 

VP() 

particle velocity V 

VPASO 

V velocity at patch 

VS2() 

temporary array 

VTRB( ) 

gas phase turbulent velocity fluctuation 

W() 

swirl velocity 

W0() 

swirl velocity of last time step 

WM( ) 

molecular weighting factor 

WP() 

particle velocity 

WTRB( ) 

gas phase turbulent velocity fluctuation 

X() 

grid point X position 

XBC() 

X projection of the first grid point away from wall 

XLEN 

X direction expansion of rectangular grid 

XP() 

particle location 

Y() 

grid point Y position 

YBC() 

Y projection of the first grid point away from wall 

YLEN 

Y direction expansion of rectangular grid 

YMAS( ) 

mass fraction 

YOP( ) 

droplet deviation from spheroidicity 

YOPDOT( ) 

time rate of change of YOP 

YP() 

particle location 

ZP() 

particle location 
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A NOVEL GAS-DROPLET NUMERICAL METHOD FOR SPRAY COMBUSTION 


C. P. Chen 
H. M. Shang 
Y. Jiang 

Department of Chemical Engineering 
University of Alabama in Huntsville 
Huntsville, AL 35899 


ABSTRACT 

This paper presents a non-iterative numerical tech- 
nique for computing time- dependent gas-droplet flows. 

The method is a fully-interacting combination of Eule- 
rian fluid and Lagrangian particle calculation. The in- 
teraction calculations between the two phases are formu- 
lated on a pressure- velocity coupling procedure based on 
the operator-splitting technique. This procedure elimi- 
nates the global iterations required in the conventional 
particle- source- in-cell (PSIC) procedure. Turbulent dis- 
persion calculations are treated by a stochastic procedure. 
Numerical calculations and comparisons with available 
experimental data, as well as efficiency assessments are 
given for some sprays typical of spray combustion applica- 
tions. 


NOMENCLATURE 

Co Drag coefficient 

Cp Model constant for the k — € model, = 0.09 

d p Droplet diameter for computational particle p 

f CdRc p / 24 

= 1 , for Re p < 1 
= 1 4- 0.15 Re° p 687 , for Re p > 1 
F{ Interaction force due to droplets 

Fpi Particle drag force 

gi Gravity vector 

k Kinetic energy 

N p Number of droplets for each computational 

particle p 

NP Number of computational particles 

P Mean pressure 

r p Droplet radius 

Rc p Particle Reynolds number 

— | 
n 

Si Source terms in the momentum equation 

t c Eddy life time 


ti nt Droplet-turbulence interaction time 

ttr Particle transit time 

U Particle relaxation time = p d d* /ISp 

Ui Instantaneous velocity for gases 

= iTi(inean) + (fluctuation) 
Instantaneous velocity for droplets 
r Effective particle relaxation time 

p Fluid (gas) density 

pd Droplet density 

p Gas viscosity 

e Dissipate rate of turbulent kinetic energy 


I. INTRODUCTION 

Numerical modeling of two- phase, turbulent react- 
ing flows has practical applications in the development 
and design of many power generating devices such as in- 
ternal combustion engines and liquid-rocket engines. In 
these flows, two-way coupling in terms of momentum, 
heat and mass transfer between underlying gas turbu- 
lence and dispersed spray droplets play one of the most 
important physical role of mixing and consequent spray 
combustion. In the last two decades computational tech- 
niques have been developed to characterize the dispersion 
of spray droplets in a gas and the influences of droplets 
on the gas dynamics. Generally, there are two approaches 
commonly used to predict gas- droplet flows. One, called 
the ” two-fluid 11 model or Eulerian approach. In this ap- 
proach, the effect of two-way coupling is incorporated as 
extra source terms in the continuum equations for both 
phases. The advantage of using a continuum approach is 
its relatively computational efficiency especially for mono- 
dispersed systems. 
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Another approach, the so called ” tracking” or La- 
grangian approach, treats the particles as discrete en- 
tities in a turbulent flow field and their trajectories are 
calculated. This approach has the flexibility of handling 
poly- dispersed spray and the two-way couplings axe usu- 
ally accomplished through the particle- source- in- cell 
(PSIC) technique with exchange of momentum, heat and 
mass between the two phases. Both approaches have been 
studied extensively [1,2] and comparative performances of 
these two approaches have also being investigated recently 
[3,4]. For typical spray combustion applications in which 
dense spray effects such as droplet collision, breakup and 
coalescence are important and drop dispersions are char- 
acterized by a non-uniform particle size distribution , the 
discrete particle approach is more convenient for repre- 
senting the poly-dispersed spray. 

In calculating turbulent gas-droplet flows, the most 
common discrete particle method is the stochastic sep- 
arated flow approach as first described by [5] utilizing 
the method of Crowe et al. [6]. In this method the liq- 
uid spray droplets are represented by a finite number of 
computational parcels and a random sampling technique 
is entailed for instantaneous gas flow properties based on 
a specified turbulence model and the resulting fluctua- 
tions are used in a Lagrangian computation for parcel 
motion. This method has been used primely for statis- 
tically steady flows in which the global iterations [7] be- 
tween the continuous phase solver and particle equation 
of motion are invoked [8,9]. For transient problems, this 
global iteration at each time step can introduce such ex- 
cessive computational requirements that numerical simu- 
lation becomes impractical. Dukowicz [10] has introduced 
a time-splitting method to couple the gas-particle interac- 
tions in a transient calculation. The numerical scheme 
used in his method is based on the SOLA code which 
utilizes a pressure substitution scheme. A similar time- 
splitting method has been used very recently by [11] in a 
pressure correction scheme for steady- state calculations 
and by [12] in a density- based scheme. 

The main effort of this paper is to present a numer- 
ical method for coupling the gas- droplet interactions us- 
ing a pressure correction scheme. This method utilizes 
the operator- splitting technique in deriving a predictor 
multi- corrector sequence which eliminates the global itera- 
tion between the two- phases at each time step. We have 
found that this method is efficient and that the required 
number of computational parcels to achieve satisfactory 
accuracy is also not excessive. In the following sections, 
formulations and validations of this method are presented. 


n. NUMERICAL MODEL 

The gas flow was formulated using the Eulerian con- 
servation equations of mass and momentum. The spray is 
described by the discrete particle method formulated on a 
Lagrangian frame. The spray is assumed to be sufficiently 
dispersed (no collision between droplets), and for simplic- 
ity, the gas flow is assumed to be close to incompressible. 


The governing equations are: Gas phase 


dp 

dt 




( 1 ) 


dpUi d 


dP 




dxj dxj 


( 2 ) 


Here all variables axe ensemble- averaged mean quantities, 
and 


r «— + + + 

in which p t is the eddy viscosity. Particle phase 
dxi 


dt 


= Vi 


(3) 


(4) 


dvi 

~ Fp* + 9i (3) 

Since the formulation here is essentially a statistical 
approach, each computational parcel represents a large 
number of droplets having equal loca tion , velocity, size 
and temperature. The two-way coupling between the two 
phases is accounted for by the interaction terms, where 


p _Ui+u'i- Vi 
Fpi ~ 

(6) 

NP 

5 m = }2N p m tv , p /dV 

p=i 

(7), 


for evaporating spray [13]; 


Fi = - |w; 3 ,N,(^) p )/dV (8), 

P ~ 1 

in which dV denotes the computational cell and the ef- 
fective relaxation time r = <*//, t* = Pdd %/ 18 and 
/ — CoRcp/24 . 

It suffices here to illustrate the calculation method 
for non-evaporating case (m etftP = 0). The method is 
based on the operator splitting technique attempting to 
reach accurate transient solution after prescribed predictor- 
corrector steps for each time-marching step. The general- 
ization of this operator-splliting technique for deriving 
pressure- correction equations suitable for all speed flows 
is described separately [14]. We focus in this paper on the 
coupling of the dispersed phase in the solution procedure. 

Discretization of the gas phase governing equation 
uses the finite volume approach. Differencing in the tem- 
poral domain employs the implicit Euler Scheme. All the 
dependent and independent variables are stored at the 
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same grid location and the variables at the finite control 
volume boundaries are interpolated between adjacent grid 
points. The discretizations have been performed on a gen- 
eral non- orthogonal curvilinear coordinate system with a 
second order upwind scheme for convection terms and the 
central differencing scheme for diffusion terms [15]. 

For incompressible gas flows considered here, the 
pressure- velocity coupling between the momentum and 
continuity equation is an important issue since density is 
constant. In the pressure-correction method, the deriva- 
tion of the pressure equation which includes the effect of 
droplets plays a key role in determining whether the ve- 
locity field satisfies the local mass continuity equation. 

In this study, a non-iterative operator-splitting algorithm 
following the spirit of [16,14] is used to derive a predictor- 
corrector sequence. We seek the finite difference form of 
the governing equations (2) and (5) as follow : 


Predictor step : 


(£ - a„)u; = H\un - a ,p n + Si - s:u: + r? (u) 


v;-v? i? + («{)"-»* 

“a T = 9i + ^ 


(15) 


The quantities 5” , R” are determined from the ex- 
isting flow fields. These values of U * and v* are used to 
evaluate r*,S * and R * such that a second approximation 
to the gas velocities can be performed : 


( At ~ Ao)U i T = ~ A,P ” + S< " S ‘ U ' T + R * (16) 

By subtracting eq.(14) from eq.(16), we also obtain the 
velocity correction equation : 


( At “ Ao)Cf * n+1 = H '( U ? +1 ) ~ A.P n+1 + Si +F^ +1 (9) 

and 

®r +1 -v? ur'+n'i-v ?* 1 , .... 

— 57“ ^ + * (10) 

The effective relaxation time scale r is evaluated 
at the second corrector level (**), to be defined later. 

The superscripts n and n + 1 denote time level t n and 
t n+1 respectively. Operator A 0 and H\ ) are constructed 
from the second-order upwind scheme for the convection 
terms and the central differencing scheme for the diffu- 
sion terms, and 5, is the source term associated with the 
Cartesian mean velocity component [14, 15]. It has been 
shown in [10] that the ensemble averaged interaction term 
F{ can be replaced by volume averaging. We split this 
term as follow : 


(| - - a„)(u; t - u n = s:u- T + s:u * +r:-r: a?) 

For the first corrector step, the momentum equation is 
approximated by 

( At “ A ° )Ur = H 'W T) ~ A<P * + S > " S * Ur + R * (1S) 
By subtracting eq.(17) from eq.(18) 

- a„ + s:wr - uf) = h\u: t - un 

-A <(P* - P") (19) 

Taking divergence of eq.(19) and invoke continuity {VU** 
= 0), the pressure correction equation is obtained : 


FT =-OT* +1 + iC (11) 

in which S ** and R ** are obtained by rearranging equa- 
tion (10) : 

- NP 

5 " = 5vE JV P m p/( A< + r D (12) 

P 


K' = ^7 E N P m p/(M + ~ + ?.At) P (13) 

P 

and m p = |7r r*pd is the particle mass. The parameter 
S? and Rl* are momentum control volume quantities 
depending on available particle information at second cor- 
rector level to be discussed later. 


A t [£>;A i (P*-P n )] = AiU: T +Ai{D:H'(U: T -U*)} (20) 


Here we have used the short notation D* u = — A 0 4- 

SJ)" 1 . The particle momentum at this level is then 


vf* - u? U* m + K) n - vf* 

— 7~ = 9i+ ^ 


( 21 ) 


The values obtained at this level are used to calculate 
r**,S** , R** and further update the velocities for gas 
phase : 


(^ - A 0 )ur T = H\u: T ) - a , p* + 5. - s: m ur T + K m 

(22) 


or, the velocity correction equation (eq.(22) - eq.(lS)) : 


By operator-splitting method, we divide the predictor- 
corrector procedure as following : 


( f j - A 0 xur T - ur ) = -s-ur T + sw 


+ 

(23) 
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At this stage, the mean velocity field satisfies the 
continuity constraint. To further satisfy the momentum 
conservation, a second corrector step is used : 


( jl _ Ao)vr = H'(ur T ) - AiP ** + Si - srur + r? 

(24) 

By subtracting (22) from (24), taking divergence and in- 
voke continuity again, we obtain the following pressure 
correction and velocity correction equations : 


A4z?r a,(p** - p*)] = a iUr T + a i \D:*H , (ur T - ^ ,* T )] 

(25) 


{ At~ Ao+ 5 ”‘ ){C7 '’*‘ - U < mT) = H ' {UrT - V ' T) 


— Aj(P** — P*) (26) 

Following [16], it can be shown that the errors introduced 
by the operator- splitting procedure is less than the trun- 
cation errors of the finite difference scheme used in the 
governing equations (9) and (10). Note that the effec- 
tive relaxation time r depends on the drag function which 
should contain the effects of turbulence. We therefore cal- 
culate (u ■)* at this stage by a stochastic method with the 
fc— € turbulence model [3,13]. A one-predictor (implicit )/one- 
correction (explicit) procedure for k and t equations sug- 
gested by Issa [11] (see also [14]) has been used in this 
study. We then let U*** and P** be the value at t n+1 
level and add the (ttj)* to update the final time level par- 
ticle velocities using the equation : 


v?+ l - v[ 

At 




+ 0 * 


(27) 


This brings all variables to the new time level. The time 
is then incremented and the new predictor-corrector re- 
peated with the new velocities. This algorithm is used 
for simulation of transient phenomena. If only a statis- 
tically steady solution is desired, then the time steps for 
gas phase and particle phase can be made unequal ; also 
the U* T and U** T calculation steps can be neglected. 


HI. NUMERICAL APPLICATIONS 

The above mentioned procedures have been coded 
into the MAST (Multiphase All- Speed Transient Flow 
Solver) program for various two-phase flow calculations. 
We first present the solid particle dispersion calculation in 
neaxly-homogenous turbulence to calibrate the stochastic 
simulation of par tide- turbulence interactions. 

Discrete Phase Turbulent Dispersion 

The stochastic technique used for modeling the dis- 
crete phase turbulent dispersion is similar to the meth- 
ods used by, for example, Dukowicz [10], Gosman and 


Ioannides [5], and Shuen et. al. [17]. The continuous 
phase turbulence is assumed to be isotropic, and the ran- 
dom turbulent velocity component u' are assumed to 
have a Gaussian probability distribution with standard 
deviation (2A:/3) 1 ^ 2 . The u[ at any required location is 
then obtained by randomly sampling the distribution 
and changes discontinuouslv after each passage of the 
droplet-turbulence interaction time tint ♦ The time t; nt 
corresponds physically to an eddy life time t e or to a 
time t tr for a particle to traverse an typical turbulent 
eddy. In [3,5,10,17] different formulations in choosing 
the t e , t tr and tint are proposed. Here we choose i t — 

1.65 C^ 4 fc 3 / 2 /e for the characteristic size of an eddy and 
t e = 3 Cpk/e for the eddy lifetime to match the asymp- 
totic dispersion analysis of Hinze [18]. 

The particle dispersion experimental setup of Snyder 
and Lumley [19] in a grid-generated turbulent flow was 
used for the numerical model validation. Particle densi- 
ties and sizes are chosen to examine the phenomena in 
which the eddy lifetime controls interaction times (46.5 
pm diameter hollow glass) , the transit time controls 
interaction times (87.0 pm corn ), or the controlling- 
interaction times undergo transition from transit time to 
eddy life time (87.0 pro solid glass). In this experiment 
fluid turbulence intensities and length scale information 
were measured. The particle calculations were started at 
the experimental particle injection point oixjm = 20 
(m is a 2.54-cm-squaxe mesh). The particle velocity was 
assumed equal to the mean fluid velocity of 6.55 m/sec. 
5,000 computational particles were sampled to calculate 
the resulting mean squared dispersion with respect to 
time. 

Comparison of the predicted and measured particle 
dispersion is shown in Fig.l. The agreement is consid- 
ered quite good. The comparison here is more favorable 
compared to the previous calculations by [17] and [5] es- 
pecially for the medium particles for which the controlling 
particle /turbulence interaction time goes through a tran- 
sition from t t to t tr - Hi this study, we did not estimated 
tint as suggested by [17] in which t tr was calculated from 
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a simplified B-B-0 equation without the gravity effect. 
Instead, we follow the stochastic procedure suggested by 
Nichols [20] and trace particle trajectories as time pro- 
gresses. This method has the flexibility of taking into ac- 
count both the gravity effect (crossing trajectory effect) 
and the non-Stokian drag law and gives more statisfac- 
torv results for medium particles. 


Although the initial jet penetration depends greatly on 
the assumed initial spray condition, the consequent good 
prediction at the latter time demonstrates the accuracy of 
the numerical scheme and importance of the droplet-gas 
interaction. 


Single Fuel Injector Case 

The experiments of Hiroyasu and Kadota [21] are 
used to validate the transient non-evaporating spray cal- 
culations. To rigorously model this flow, detailed atom- 
ization processes have to be resolved at the injector noz- 
zle. Since this phenomena is not modeled in this study, 
information about the injected droplet size distribution, 
as well as the velocity distribution has to be estimated 
in vicinity of the injector at which the calculation starts. 
Following the suggestions of [10], the particle initial in- 
jecting velocity was determined by the mass flow rate and 
pressure difference of the nozzle. The estimated initial 
spray angle was also guided by [10]. The initial particle 
size distribution was given by the following form : 

fi{d p ) = -jr-exp (-^) 

Uzi U 32 

where D 32 Is the Sauter mean diameter. 

The test conditions Eire given in Table 1. The spray 
was assumed to be axially symmetric, and the calcula- 
tion was carried out in cylindrical coordinates. The com- 
puted penetration of the tip of the spray as a function of 
time for the two test conditions are qualitatively shown 
in Fig. 2 and quantitatively compared with the experi- 
mental data in Fig. 3. The comparison is very good here. 


FM 


>02 ms 


P=30 atm. 


*M)5 ms 1.0 ms 

Figure 2. Particle Plots of a Single Orifice Spray. 



TABLE 1 

SINGLE-ORIFICE INJECTION PARAMETERS [21] 


Chamber 

Injection 

Gas 

Mass 

Sauter Mean 

Nozzle Pressure 

Gas Pressure Velocity 

Density 

Flow 

Radius(SMR) Difference 

(atm) 

(m/sec) 

{kg/m 3 ) 

(kg/sec) 

pm 

AP(atm) 

1 

122:2 

1.123 

0.00726 

5.0 

98 

30 

102.5 

33.70 

0.00609 

5.0 
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Fuel: Diesel fuel oil, pd = 840 kg/m 3 
Ambient Gas: Nitrogen 
Nozzle Diameter: 0.3 mm 
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PENETRATION 



Figure 3. Comparison of Computed Spray Penetration 
with Experimental Data. 


Hollow-Cone Spray 

A polydispersed pulsed hollow-cone spray case of 
practical importance is also chosen for the test condition 
listed in Table 2. Fig. 4 shows the particle plot and the 
gas velocity vectors for a 30° spray. With the back pres- 
sure 1 atm, the interaction between the gas and droplets 
is seen to be rather strong. The shape of the spray is 


no longer conical even for very short time and the spray 
penetration is suppressed due to the interaction of the 
droplets with the induced air flow. The gas velocity vec- 
tors indicate the presence of a vortex near the head of 
the spray, which curls the spray tip toward the outside of 
spray. A substantial region of strong inward flow in the 
center of the cone near the injector was also observed. 
These flow patterns and spray shapes compared quite fa- 
vorably with the experimental observations (c.f.[22]). 

The efficiency assessment of the present numerical 
method was shown in Table 3 for the single-orifice and 
hollow-cone spray cases. The CPU times on a CRAY 
X/MP using the MAST code utilize the preset method 

and the TEACH code with PSIC method [3] for both 
transient spray calculations using At = 0.1ms are given. 

It can be seen that the amount of CPU- time is reduced 
about one order of magnitude using the present calcula- 
tion procedure. Also, the present method is rather par- 
ticle number independent. This is due to the fact that 
particles are injected at each time step and the source 
terms in the continuous phase axe updated for all the 
particles at each Eulerian control volume. While in the 
TEACH /PSIC method, all the particles have to be tracked 
and the continuous phase flow field is held frozen between 
the global iterations at each time step. This PSIC algo- 
rithm thus requires substantial computer time and is in- 
herently unsuitable for transient calculations. 


TABLE 2 

HOLLOW-CONE SPRAY PARAMETER 


Chamber Injection 

Injection 

Gas 

Mass 

Sauter Mean 

Gas Pressure Velocity 

Angle 

Density 

Flow 

Radius(SMR) 

(atm) (m/sec) 

(Degree) 

(kg/m 3 ) 

(kg/ sec) 

fim 

1 20.0 

30 

1.123 

4 x 10~ 4 

2.5 


TABLE 3 



EFFICIENCY ASSESSMENT (CPU Time) 



MAST-2D 

TEACH/PSIC 

SINGLE-ORIFICE SPRAY Particle # 


Particle # 


41 x 61 Grid 

600 

126.9 sec 

800 

1420 sec 

300 Time Steps 

1200 

135.7 sec 



HOLLOW-CONE SPRAY 

400 

74.9 sec 



31 x 31 Grids 



800 

934 sec 

200 Time Steps 

1000 

88.3 sec 
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Figure 4. Velocity Vectors and Particle Plot of a 30° Hollow-cone Spray. 


IV. SUMMARY 

A new numerical scheme based on a non-iterative 
predictor-corrector pressure velocity coupling procedure 
has been developed for transient gas-droplet two-phase 
flows. The present scheme is formulated on a Euler ian - 
Lagrangian analysis and the two-way interaction between 
the two-phase is handled through a strong coupling pro- 
cedure. This procedure eliminates global iterations con- 
ventionally used in the PSIC procedure and shows dras- 
tic saving in CPU time for transient spray calculations. 


This method has been extended to evaporating spray 
calculation recently [13]. Good agreements between the 
calculated results and the experimental data have been 
obtained despite the uncertainties in the inlet conditions 
for both fluid and droplets. The development of the tech- 
nique is based on the assumption of dilute ( non-iterating 
) sprays. The computer code, however, allows easy alter- 
ation of models, so that an appropriate model to suit the 
physical problem of interest can be quickly implemented. 
Experimental studies with better defined inlet conditions 
would be extremely useful in further model validations. 
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